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COMBUSTION OF HYDROGEN-AIR JETS 
• IN LOCAL CHEMICAL EQUILIBRIUM 
(A Guide to the CHARNAL Computer Program) 


by 

D. B. Spalding, B. E. Launder, A. P. Morse and G. Maples 


of 

Fluid Mechanics and Thermal Systems, Inc. 
Route 2, Box 11 
Waverly, Alabama 36879 


1. Introduction 

The computer program CHARNAL (Calculator of Hydrogen-Air Reactions 
for NASA Langley) generates finite-difference predictions of turbulent, 
coaxial hydrogen-air jets undergoing combustion. The jets may be free, 
as indicated in Figure la (in which case the external stream is assumed 
to extend to arbitrarily large radius from the symmetry axis) or enclosed 
as in Figure lb. At any point in the flow t;he mass fraction of the con- 
stituents of combustion is found on the assumption that chemical equili- 
brium prevails, the constituents being H 2 , 0 2 , H 2 0, 0, H, OH and N 2 . 

The present report details the mathematical and physical basis of 
CHARNAL, discusses some sample predictions and provides a guide to the 
computer program itself. Section 2 is concerned with the first of these 
aspects: the basic conservation equations of momentum, stagnation 

enthalpy and chemical species are presented first and these are followed 
by a description of the turbulence and combustion models employed. A 
user's guide to the computer program appears in Section 3 while defini- 
tions of FORTRAN symbols and a listing of the program itself are contained 
in the Appendices. Thereafter, Section 4 presents and discusses the out- 
come of some test cases and, finally, Section 5 suggests some directions 
that further developments to the CHARNAL program might take. 



2. The Mathematical and Physical Model 


2.1 The Mean Flow Conservation Equations 

CHARNAL calculates the steady state distributions within the jet of 
mean streamwise velocity, temperature and mass fraction of elemental hydro 
gen by reference to the conservation lav/s of momentum, energy and chemical 
species. These laws are expressed in terms of the following set of para- 
bolic partial differential equations expressing respectively the transport 
of streamwise momentum, stagnation enthalpy and hydrogen mass fraction. 









( 2 . 1 - 1 ) 


pu + PV 111 = 1 1- 

3x 3r r 


P t r + 1 1_ 
t dr r 3r 


( r k- r h) 


~ + z. I r .- r h 
dr J \ o h 


K- r h) 


ar 


(u 2 /2) 


K 


3m. 

_2 

3r 


( 2 . 1 - 2 ) 


pu — + pv — = — — r r — 
3x 8r r 3r m 3r 


(2.1-3) 


The mass fraction of elemental hydrogen, f, is used as a variable for the 
reason that, unlike H ? , even during chemical reaction it remains a con- 
served property. The method of determining the individual chemical 
constituents of the flow is described in Section 2.3. 


The set of partial differential equations is completed by the continu- 
ity equation in which the streamwise and radial velocities are connected 
by: 
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In fact CHARNAL solves the parabolic transport equations fcast in a 
Von Mises system of coordinates (i.e., x and stream function as independent 
variables). Tiiis transformation has the effect of eliminating the radial 
velocity v from the equations, and,' hence no explicit recourse needs to 
be made to equation (2.1-4). 

*>' The temperature of the mixture, T, is obtained from known values of 
h, u and the mass fractions of the chemical constituents of the mixture 
from the expression: 



(2.1-5) 
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2.2 The Turbulence Model 



and applied to numerous turbulent free shear flows in Reference [2]. 
According to. this model the magnitude of the viscosity depends only on the 
local values of the turbulence kinetic energy, k the dissipation rate of 
turbulence energy, e and the fluid density. They are connected by the 
formula: 


= C u pk 2 /e 


( 2 . 2 - 1 ) 


The quantities k and e are found by way of the following pair of transport 
equations which are both similar to (and solved simultaneously with) those 
governing the mean flow: 
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The quantities C , C-j , and 0 E are dimensionless and are given 

the constant values belcw 
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s. 
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These values are the same as those recommended in [2]. Equations (2.2-2) 
and (2.2-3) do not provide a physically exact prescription for finding 
k and e. Such a prescription is not possible because the exact equations 
for k and e contain correlations whose magnitude is not directly ascer- 
tainable. These correlations are therefore approximated in terms of k,e 
and the mean velocity field; the approximated terms in (2.2-2) and (2.2-3) 
are the ones with the empirically determined coefficients. Models of the 
above kind, while being sufficiently simple not to affect computer costs 
significantly, have been found (see for example, references [2] and [3] ) 
to possess considerable width of applicability, precisely the same model 
predicting features of both wall and free turbulence. It is probably the 
best model available at present for the kind of shear flows that CHARNAL is 
designed to compute. 

The transport coefficients in the hydrogen-element and stagnation- 
enthalpy equations are given by 


T h'~n/% '» r m = H /a m 


(2.2-4) 


In the free jets and in the confined jets (provided the jet has not spread 
to the pipe wall) 


o h = o ;n = 0.7 (2.2-5) 

Once the jet has filled the pipe the effective Prandtl/Schmidt number is 
obtained from the formula: 


o h = o m = 0.95 -.45 (y/R) 2 (2.2-6) 
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where y is distance in the radial direction measured from the pipe wall. 
The above variation, proposed by Rotta [4] and used by several workers 
since, is generally in accord with experimental data of the turbulent 
Prandtl number in fully-developed pipe flow. Note that because the same 
numerical values are assigned to and r m , one of the source terms in 
equation (2.1-2) vanishes. 

2.3 The Combustion Model 


The equilibrium composition of the hydrogen air mixture can be cal 
culated by reference to the set of reversible reactions: 


H 2 z=± H + H 

(2.3-1) 

o 2 ^ 0 + 0 

(2.3-2) 

OH ^ 0 + H 

(2.3-3) 

H;,0^ OH + H 

(2.3-4) 


The relative mass fractions of the above constituents are found by presum- 
ing that chemical equilibrium prevails at each point in the flow. Thus: 


m^/m^ = (2.3-5) 

m °2^ m 0 = (2.3-6) 

m Ol/ m O m H = 1 '3 (2.3-7) 

p H 2 0 /m H m 0H - K 4 (2 - 3 - 8) 


where the K's are functions of temperature and pressure. In addition we 
have by definition 


and 
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In the above equation set, the quantities f and X are to be considered 
kncwri: f is found from the conservation equation for total hydrogen 

(equation 2.1-3) and X is simply OFAC (1-f) where OFAC is the mass .frac- 
tion of oxygen in the external stream. The equilibrium 'constants 
K-| K* are functions of temperature and pressure. CHARNAL incor- 

porates tne dependencies proposed by McBride [5]. 


Equations (2.3-5) - (2.3-10) thus provide a set of six equations in 
the six unknowns m.. , m u , m.. m A1J , m A , m A . The equations are non-linear 

n HgU Un v» 

(indeed, highly non linear for regions of flow where the fuel/air ratio 
is nearly stoichiometric); their solution must therefore proceed itera- 
tively. The solution technique adopted is described in Section 3-4. 


is s 


The mass fraction of N« (which is considered to be entirely inert) 
imply (1-X-f). 


Because the flow is turbulent, the level of hydrogen and of other 
variables will be continuously fluctuating about their mean values. The 
magnitude of the mean square hydrogen fluctuations, g, is found from the 
following transport equation developed and tested by Spalding [6] 
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where the constants o g , C^ , and are assigned the values of 0.7, 2.8 

and 2.0 respectively. Solution of the g equation enables the region of the 
flame over which stoichiometric condi tons occur to be calculated. The 
inner boundary to this region occurs at a radius from the axis where f-g^ = 
^stoich and outer b° unc * ar -7 Vvhere f+gh - fgtoich* The predicted path 

of these two surfaces for a hydrogen jet burning in an air stream is shown 
in Figure 10. 

3 . 0 Details of the CHARNAL Computer Program 

CHARNAL is a custom-developed version of the Patankar-Spalding PASSA 
program which, in turn, is a newer, more economical and flexible version 
of GENI1IX published in Reference [1]. 

Section 3.1 below gives a brief summary of the general Patankar and 
Spalding method of solving the boundary layer equations. Section 3.2 
lists the major differences between CHARNAL and the original GENMIX pro- 
gram, while Section 3.3 discusses in detail the listing of the present 
program (which is to be found in Appendix 2). Section 3.4 is devoted to 
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a discussion of subroutine LAi'IGLEY , which has been developed exclusively 
for the current work. Finally Section 3.4 gives overall instructions for 
the running of the program, in particular the input information that the 
user must provide. 

3.1 The Patankar-Spaldlng Boundary Layer Procedure 

The main features of the computational method are as follows: 

(i) The primary differential equations (i .e. , the transport equations for 

A . . ’ • ... 

u, h, f, g, k and e) are transformed so that the independent variables 
are the longitudinal distance, x and the dimensionless stream function, 
u) defined as: 


u = (f - ’• 'Pf) (3.1-1) 


where ij; itself is obtained from the relation. 


3t!/ 



(3.1-2) 


Subscripts I and E refer respectively to the "inner" and “cuter" 
boundaries of the flow; for axisymmetric flows the inner surface is 
always the one nearer the symmetry axis. The resulting differential 
equations are all of the basic form, 
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where the terms on the left-hand side represent convection by the mean 
flow and those on the right-hand-side express respectively the dif- 
fusion and source of the entity The coefficients, (a) and (b) 
are functions of the entrainment rates, whilst (c) involves the ef- 
fective diffusion coefficient. 


(ii) The differential equations are expressed as finite-difference equa- 
tions, connecting the values of the dependent variables which prevail 
at the intersection points of a grid defined by lines of constant x 
and to. These finite -difference equations are formed by integrating 
the differential equations over small control volumes associated with 
each grid node. 
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(iii) The integration proceeds by "marching" downstream, the values of the 
variables at grid points at the next downstream station being calcu- 
lated from those at the upstream station and placed into the array 
locations occupied by the latter. At each forward step, new values 
are ascribed to •!.•£ and tpj, the stream functions at the grid boundaries. 

These values, together with the continuity equation, determine the 
geometrical location of the boundaries; and this determination is so 
arranged that the boundaries enclose all the fluid having a signifi- 
cant value of the shear stress (or other flux) without enclosing 
appreciably more than this. This feature allows the method to achieve 
good numerical accuracy without employing an excessively fine grid. 

( i v ) The finite-difference equations are formulated in an implicit manner, 
and solved by means of the well-known algorithm for tri -diagonal 
matrices. This allows large forward steps to be made without insta- 
bility. The equations are linearized, upstream values of the trans- 
port properties being supposed to prevail over the whole of the 
forward step. 

(v) The source terms are usually (but not necessarily) expressed as linear 
functions of the upstream and downstream values of the dependent vari- 
ables. For example in the turbulent kinetic equation. 


Source of (k) = 

(§? - 

Hpe/k) 

k D 

(3.1-4) 
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U 
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where the subscripts IJ and D. denote upstream and downstream values. 
This practice allows large forward step sizes to be used without the 
onset of inaccuracy or instability. 

For further information on the general structure of the solution pro- 
cedure the reader is referred to Reference [1]. 

3 . 2 New and Improved Featur e s of the PASSA family of programs 

PASSA programs retain the general width of applicability of the basic 
procedure, but are arranged to be more economical in terms of both execu- 
tion time and storage capacity. The principal differences between the new 
method and its predecessor may be summarized: 

(i) In PASSA, the finite-difference equations are solved sequenti ally . 
That is, for any variable, the coefficients of finite-difference 
equations are formed and the equations solved before moving on to 
the next differential equation. This contrasts with the former 
procedure where the coefficient for all the differential equations 
were stored simultaneously and then the matrix was inverted for each 
equation in turn. 
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(ii) The PASSA program does not make use of "slip nodes" in the formulation 
of the difference equations at the boundaries. In the earlier program, 
these nodes were employed to obtain the correct gradients of the depen- 
dent variables at the edge regions. PASSA treats the nodes near each 
flow boundary in exactly the same way as the interior nodes. 

(iii) In PASSA all di ffusi vities and gradients of the dependent variables, 
such as the viscosity and the velocity gradient are evaluated at the 
actual grid nodes, whereas in the earlier program these quantities 
were evaluated at the edges of the control volumes (which by definition 
were midway between individual nodes in terms of w) . PASSA also con- 
tains the option of using castellated profiles instead of linear vari- 
ations of <p between the grid nodes. 

3 . 3 Details of the CHARNAL Program 

A flow diagram for CHARNAL is provided in Table 1. The program con- 
sists of a MAIM program and a number of subroutines of which the most 
important are AUX, STRIDE, OUTPUT and LANGLEY. 

3.3-1 The MAIN Program 

MAIN contains the starting and stopping points of the computation and 
communicates directly or indirectly with all the other subroutines. It 
comprises twelve "chapters", each performing a specific function in the 
computational procedure. The most important of the operations are men- 
tioned below: 

Chapter 1. Values are assigned to various indices which control aspects 

of the computation throughout the program. The main categories 
are: 

(i) the specification of the number of grid nodes and 
differential equations to be solved. 

(ii) the nature of the flow boundaries (wall, free boundary 
or axis of symmetry) 

(iii) the control of input, and 

(iv) the designation of the flow type (number of chemical 
species, nature of initial conditions, etc.). This 
topic is discussed in more detail in Section 3.4. 

Chapter 2. This chapter selects the primary dependent variables and 
auxilary quantities to be calculated in the program. Com- 
ment cards have been included here for the user's benefit. 

Chapter 3. Material constants such as molecular weights and the universal 
gas constant, turbulence parameters and Prandtl/Schmidt 
numbers are assigned values. The S. I. system, of unit is 
used throughout the program. Subroutine LANGL1 
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is called to provide data for the enthaplies of the chemical 
species and the equilibrium reaction constants. 

Chapter 4. Specified here are the flow geometry (plane or axi symmetric) , 
inclination of the streamlines to the axis of symmetry and 
the cross-stream distances for the initial profiles. 

Chapter 5. The initial profiles are read in from data cards. Values 

are assigned to the axial velocity, absolute temperature and 
species concentrations at the grid nodes. ATI the data input 
are in dimensionless form, having been normalized with the 
largest value of each dependent variable at the initial 
station. In developing CHARMAL, fifteen test cases (speci- 
fied by NASA Langley) have been run v/hich required two dif- 
ferent types of initial-profile specification. 

( i ) Cont i nuous profiles (Test Case nos. 1-10,13) 

Figure 2a shows the initial-velocity-prof ile typical of 
these test cases. The two streams are separated by a 
wake region caused by the interaction of the boundary 
layers on the dividing wall. Since the velocity is 
uniform near the jet centreline, the computation starts 
from a mixing-layer region with entrainment at the 
inner flow boundary until this boundary grows to the 
axis of the jet (i.e., the end of the potential core). 

( i i ) Step profiles (Test Case nos. 11-12, 14-15) 

In these four test cases, the boundary layers on the 
dividing wall are ignored and a step-change in the 
velocity profile is assumed (see Figure 2b). In these 
cases, the computation starts from a very thin mixing 
layer (with an assumed linear velocity profile) in the 
immediate vicinity of the step-change. 

The initial profiles of turbulent kinetic energy are 
then evaluated by assuming a constant ratio with the 
shear stress (as expressed by the mixing length hypo- 
thesis), whence 



(3.3-1) 


The dissipation rates are then given by the Prandtl- 
Kolmogorov relationship 


. 3/2 . 

e = k /£ 


(3.3-2) 
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Both the mixing length and dissipation length scale 
are assumed to be proportional to a typical width of 
the shear region. The scheme adopted is shown in 
Figure 2. 

The free-stream turbulent kinetic energy is taken as 

4 x 10 of the square of the free-stream velocity. 

The same constant is used to determine the initial pro- 
file of concentration fluctuations from the local con- 
centrations of hydrogen element, i.e., g/f 2 = 4 x io~^. 

A variable of major importance in the computational 
procedure is OFAC, which is the ratio of oxygen element 
to nitrogen in the outer stream, since this ratio is 
assumed to be constant across the flow. 

Chapter 5 also calls subroutine LANGL2 to evaluate the 
initial enthalpy profile. 

Chapter 6 . The dimensionless stream function array is filled by inte- 
gration of the profiles of density and velocity according 
to equations (3.1-1) and (3.1-2). The density is obtainable 
from the ideal gas relationship, 


P = plJ/RT (3.3-3) 

where the mean molecular weight of the mixture, W is given 
by, / 



Subroutine STRID1 is called to evaluate useful quantities 
relating the individual W's. 

Chapter 7. This marks the starting point of the main computation; it 
is the point to which control is returned after the execu- 
tion of each forward step. The most important functions 
of this chapter are, 

(i) to call LANGL3 , which employs the chemical reaction 
constants relevant to the upstream conditions, (i.e., 
pressure and temperature) to determine the local mass 
fractions, 

(ii) to call LANGL4 to calculate the temperatures corres- 
ponding to these concentrations (by way of the up- 
stream specific heats) and • 
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Chapter 8. 


Chapter 9. 


Chapter 10. 


Chapter 11. 


(iii) to evaluate the density profile from the local temp- 
eratures and mass concentrations 

This chapter performs two main tasks. First it fixes the size 
the forward step and secondly it calls subroutine STRID2 to 
calculate transverse distances. The forward step is usually 
made proportional to the width of the flow, with the constant 
of proportionality small in the initial region to avoid in- 
stabilities at the start of the calculation. For confined 
flows, the streamwise pressure gradient, which is a source 
term in the axial momentum equation, is not known a priori . 
CHARNAL adopts the same non-iterative practice employed in 
GErltUX, (see reference [1]). An estimate of the pressure 
change to be experienced over a forward step is obtained by 
reference to a 1-dimensional analysis. This usually results 
ir. the area of the flow differing from the pipe cress 
sectional area at the end of the forward step. However by 
adjusting the level of dp/dx ov er the next step the differ- 
ence in area can be kept negligible (typically 0.01% of the 
pipe area). 

This chapter fixes the conditions at the flow boundaries. 

Only when a boundary is a wall- must information (either the 
value of if> or its diffusional flux) be specified at this 
point. (In .the case of a free boundary, the relevant infor- 
mation is provided in STRID3 (based on the free-stream source 
terms) while at an axis of symmetry, the zero gradient con- 
dition usually applies). 

The first chapter of AUX is called to determine the effective 
viscosity (regarded as the sum of the turbulent and I;.minar 
viscosities) at each node and to formulate the source term 
based on the axial pressure gradient. For a free boundary 
the entrainment rate is calculated via the degenerate form 
of the conservation equation for whichever of the dependent 
variables shows the largest changes near the edge of the flow. 
The entrainment is subject to certain controls to prevent the 
formation of 'tails' to the profiles and to prevent the onset 
of instability. 

This chapter deals with the output of information. At the 
first step, subroutine OUTP 1 is called to print-out inform- 
ation regarding the initial conditions. At certain axial 
positions (designated in Chapter 1) OUTP 2 is called to 
print-out such quantities as the entrainment rates, jet 
spreading rate, centreline values of the velocity, tempera- 
ture and species concentrations and the fluxes of the dep- 
endent variables. Also OUTP 3 may be called (not necessarily 
at the same stations) to print-out the profiles of quantities 
of interest (velocity, temperature, concentrations, etc.). 
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Chapter i2. The last chapter terminates the execution after a specific 
axial distance has been covered. Otherwise, control is 
returned to Chapter 7. 

3.3- 2 AUX, STRIDE and UF 

Only brief description is provided of these general-purpose subrou- 
tines as they are similar in structure and function to the correspondingly 
named subroutines in GENMIX (reference [1]). Subroutine AUX is called 
initially from Chapter 10 of r'AIN to provide the effective viscosities at 
the grid nodes. Chapter 2 is subdivided into five parts, one for each of 
the dependent variables other than the velocity. Each section evaluates 
the appropriate effective diffusivity and formulates the source arrays and 
is called in turn from STRIf 3 as the finite difference equations are 
solved sequentially. The source terms for turbulent kinetic energy, dis- 
sipation rate and mean square concentration fluctuations are linearized 
according to equation (3.1-4), while that for stagnation enthalpy is 
loaded entirely into the upstream array. 

Subroutine STRIDE is divided into four parts, of which the first two 
are largely preparatory, while the latter two contain the core of the 
numerical method of the program. STRID 0 merely sets to zero arrays 
such as those for the dependent variables and auxilary quantities. It 
is here that the decision is made as to whether to employ linear 
(KAST (J) = 0) or castellated (KAST (J) = 1) profiles between the grid 
nodes. STRID 1 evaluates useful relationships between the values of w 
at neighboring nodes, these of course remaining constant for the whole 
of the calculation as the grid is always constrained to lie between limits 
of 03=0 and w=l. STRID2 calculates the cross-stream distances at each 
axial station with but minor differences from the original program [1]. 
STRID3 contains the basic finite-difference formulation and technique 
of solution of the differential equations. Although different in appear- 
ance from the version published in reference [1], the differences are 
largely only ones of arrangement. It is not proposed to dwell on points 
of detail here, but suffice it to say that STRIDE is a subroutine which 
the user has very seldom any need to change. STRID3 terminates by 
determining conditions at the flow boundaries and by initiating the 
forward step. 

Subroutine IJF provides wall -function relationships to relate the 
fluxes (diffusional and convective) through the wall with the values of 
the dependent variables at the near-wall nodes. As appropriate to the 
high-Reynolds-number form of the turbulence model employed in this work, 
the wall functions are based on the assumption of a log-law velocity 
profile in the fully turbulent region of the flow. 

3.3- 3 OUTPUT and related Subroutines 


Subroutine OUTPUT is divided into three parts, concerned respec- 
tively with the print-out of initial values of interest and station and 
profile variables. It communicates with PROFIl and thence with subroutine 
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PLOTS, which together may be used to provide non-dimensional plots of the 
profiles. 0UTP1 reads alphanumeric information from data cards to provide 
the headings for the print-out. The meaning of the quantities printed 
from 0UTP2 and CUTP3 is given in Appendix 1. Subroutine QUTP2 performs 
the useful function of checking for the overall conservation of fluxes of 
the dependent variables; in this analysis the concentration of hydrogen 
element is a conserved property as are the axial momentum and stagnation 
enthalpy (but for free flows only). 

There is a large number of comment cards in OUTPUT, PROFIL and PLOTS 
to assist the reader in understanding the interlinkage between these sub- 
routines. PROFIL is used to normalize the profiles and communicates with 
PLOTS which scales both the abscissa (transverse distance) and ordinate 
U values) of the dimensionless plots into the range 0 to 1, (with nega- 
tive values printed as zero). Either the normalized (IPROF = 1) or the 
full dimensional (IPROF = 2) values of the dependent and auxiliary vari- 
ables may be printed out in 0UTP3. If the former approach is adopted, 
the full dimensional values are still printed at the first and last nodes 
(designated 1 and UP! respectively). The cross-stream distances Y(I) are 
treated In this way, the first quantity printed being the radius of the 
internal flew boundary and the last the distance between the internal 
and external boundaries. 

Subroutine YINT may be called at any point in the program and is 
merely an interpolation subroutine. It is useful in the determination 
of such quantities as the half-width of the jet or locating the exact 
position where the value of an entity <p is a linear combination of the 
values of 4 > at the inner and outer boundaries of the flow. 

3.3-4 Subroutine LANGLEY 


Subroutine LANGLEY has been developed exclusively for the present 
work and is subdivided into four parts, whose purposes are respectively 

(i) the loading of the individual species and the chemical 
equilibrium constants for the stipulated reactions. 

(ii) the evaluation of the initial enthalpies by simple 
interpolation amongst the input data. . 

(iii) the calculation of the species mass fractions by solu- 
tion of six simultaneous equations. 

(iv) the evaluation of the cross -stream temperatures from 
the enthalpies and the upstream specific heats. 

Data for enthalpies and chemical-reaction constants are those pro- 
posed by McBride [5]. The enthalpy of each species may be written in 
the form 


14 



(3.3-5) 


h ' h ref " / 


C dT + Ah., T f 
p f ref 


ref 

where Ahf, T denotes the heat of formation of the species at the refer- 
ref 

ence temperature. It is usual to take the datum of enthalpy as zero at 
a reference temperature of 298. 1 5° K, whence 


h = 


I 


298. 15°K 


CpdT + Ahf ’ 298.15°K 


(3.3-6) 


/ V 

0°K 


‘ hf, 298.15”K ' (h 298.15"K" h 0”K 1 


.J 


= C p T + H q (3.3-7) 

T 

where C (= 1/ T J C p dT) denotes a mean specific heat and H q is the 

0°K 

composite of the remaining terms. Values of h are read in at intervals 
of 100°K for temperatures of up to G000°K, together with the value of H 

for each species. The mean specific heats are then evaluated from equa- 
tion (3.3-7) and stored (see Appendix 1 for array locations). 

LANGL2 is called to obtain the initial enthalpy profile by linear 
interpolation of the h/T data. The enthalpy of the mixture is then 
given by 


h = 7 , m^h - 
o J 

1 m J C P,j T + 1 (3.3-8) 
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LAMGL1 is also usee! to score data for the chemical equilibrium for 
the four reactions involved in the combustion process. These are again 
taken from McBride [5] and are in the form of constants relating the 
partial pressures of the individual species, i.e., for the reaction 


aX + bY = Z 


(implying the formation of, 1 unit of substance T from (a) units of X 
and (b) units of Y), the partial pressures of the species are linked by 
a constant, K given by the relation 

p 


K p • P Z /P x a P v b (3.3-10) 

[Because K (which is a function of temperature) varies markedly, numbers 
on the data cards refer to 1 og-j Q Kp]. Kp's ma y be transformed into con- 
stants relating the species mass fractions by use of the ideal gas law, 
whence 

K m i ,n 2 /m x a m y b 

* K p (pW) a+b 'V (wh'/j (3.3-11) 

where V!. denotes the individual molecular weights and W is the mean mole- 

v . •« L 

cular weight of the mixture. Values of (KpW z /W x “Wy U ) are 

stored for each reaction involved. In using Equation (3.3-11), it is to 
be remembered that the total pressure, p must always be expressed in 
atmospheres . 

LANGL3 evaluates the concentrations of the individual species by 
six simultaneous algebraic equations, viz. 


; = m n /m 2 

m, 1 V m H 

(3.3-12) 

in, 2 = V"o 

(3.3-13) 

in, 3 ° ”0H /m 0 m H 

(3.3-14) 
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(3.3-15) 


^,4 “ m H 2 G /m H I ’oH 


H 0 W 

x = % + m 0 + “ m H 2 0 + — m 0H 


H 2 0 


OH 


(3.3-16) 


W 

"u ».i 

n 2 "H 

f - m + m u + — m + m nu 

H H w H 2 0 u OH 

2 H_0 1 OH 


(3.3-17) 


in which the equilibrium constants (which are functions of pressure, temper- 
ature and concentrations) are based on the upstream conditions, since 
the temperatures can only be calculated once the mass fractions have been 
found. 


The first function performed by LAMGL3 is to evaluate the reaction 
constants and mean specific heats by linear interpolation amongst the 
input data. Equations (3.3-12) to (3.3-17) are then solved iteratively 
until convergence to within 1/6 is obtained for all species present to 
an extent of more than 10'^, or until the number of iterations exceeds 
a certain limit. The actual iterative scheme is of necessity rather 
elaborate since the relative concentrations of the species can vary 
markedly. The following practices are adopted 

(i) 7 he concentration of nitrogen is easily calculable 

as (1-OrACMl-f). 

(ii) The concentrations of the remaining species ati the 
previous axial position are used as a 1 first estimate 
of the downstream yalues. 

1 , 

(iii) The equations for f and X[= OFAC(l-f)] are examined to 
determine the largest term on the right-hand side of 


each. For example, suppose that 





and 


nin are the largest quantities. 
2 
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( i v ) Using the values of m^ q and , equations (3.3-12) 

to (3.3-15) can then be solved to determine the con- 
centrations of the remaining four species. 

(v) m^ q and m^ can then be re-evaluated from Equations 

(3.3-16) and (3.3-17). 


(vi) The right-hand sides of the X and f equations are then 
re-examined and the next iteration cycle made to operate 
on the two largest quantities. 


(vii) Iteration is repeated until the convergence criteria 
are satisfied or the number of iterations becomes 
excessive. 


Convergence is usually rapid (2 iterations), but instability occurs 
whenever two of the terms in either of the additive equations [(3.3-16) 
and (3.3-17)] become of approximately equal magnitude. This happens 
(a) at high temperatures (> 3000°K) where the concentration of atoms 
approaches that of the molecules and (b) in the region of the stoichio- 
metric point where the concentration of combustion products becomes-large. 
In the flows under investigation, the maximum predicted temperature was 
approximately 2600°K; (b) was therefore the main source of instability. 

It often happened that on successive iterations different constituents 
were the largest terms in Equations (3.3-16) and (3.3-17). Oscillatory 
rather than convergent behavior would then result. However, under- 
relaxation of the species concentrations between successive iteration 
cycles generally achieved the required degree of convergence within 12 
iterations. A warning message is printed out whenever the iteration 
process fails to converge within this limit. 

When the species concentrations have been found, LANGL4 is called 
to evaluate the cross-stream temperatures. These are obtained from 
enthalpy profile via the relationship. 


T = (hr - Em.:H n .)/m.C , (3.3-18) 

J 0»J J Pj 

for which the mean specific heats are evaluated in LANGL3 on the basis 
of the upstream temperatures . One could regard this value of temperature 
as a first estimate of the actual value, since a more accurate estimate 
could be obtained by re-evaluating the K^’s to obtain new values of the 
constituent mass fractions, a new mean specific heat, and from (3.3-18) 
a new T. LANGLEY does not incorporate such an iteration cycle (though it 
would not be difficult for a user to add). Instead, under-relaxation of 
the temperature between upstream and downstream stations damps out any 


18 



temperature spikes which would otherwise appear at positions close to 
the stoichometric point. 

3 . 4 Use of the Present Program 

CHARNAL solves 6 simultaneous partial differential equations and 
employs typically 25 nodes (although the storage blocks have been dimen- 
sioned to allow for as many as 40). The program is written in basic CDC 
Fortran language and requires approximately 12 seconds compilation time 
On a CDC 6600 machine (FUN compiler). Linear profiles are assumed for 
the variation of the dependent variables between the grid nodes and a 
typical forward step size of 0.1 times the local shear layer width is 
employed. Approximately 15 axial steps can be executed per second. Con- 
servation of the individual fluxes is generally good to within + 0.1%. 

In the present work, the profiles have been printed out in full dimen- 
sional form and are plotted at only one axial position (corresponding to 
the end of the potential core region). The program has a modest storage 
requirement, needing only approximately 24000 decimal storage locations. 

The remainder of this section is concerned with the information that 
the user has to provide via data cards in subroutines MAIN and LANGLEY. 

In Chapter 1 of MAIN, control indices which have to be set are: 

( i ) Details of the grid and nature of the flow 

KASEN0 : Test case number 

NEQ : number of differential equations 

M : number =ef gri d nodes 

KASE i type of flow; = 1 (free jet), = 2 (confined 
flow) 

KIN : specification of internal flow boundary) 

= 1 (wall) 

= 2 (free) 

KEX : specification of external flo w boundary) 

= 3 (axis of symmetry) 

KONFIN : denotes presence of confining duct wall; 

= (1 wall present), = 2 (no wall) 

IN2, 102, IH20 : denote presence of individual species; 

= 0 (not present), = 1 (present) 

NR : number of reactions 

NS : number of chemical species 
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INUF : nature of initial profiles = 0 (continuous) 

= T (step-change) 

: initial axial distance (in nozzle diameters) 

XULAST : last axial distance (in nozzle diameters) 

YIN : radius of internal flow boundary (in nozzle 
diameters) 

YOUT diameter of jet nozzle 

RDUCf : radius of confining duct 

PRESS : static pressure (in atmospheres) 

■ • TA : maximum temperature at initial station 

j'UIH : maximum velocity at initial station 

YW1 : extent of inner shear, region (in nozzle 
diameters) i.e., distance -from axis of 
symmetry to minimum in velocity profile 

UWl : velocity at position denoted by YIJ1 

IU1 : grid node corresponding to YW1 and UHl 

TDUCT : temperature of duct wall (for confined 
flows only) 

1 ' ( i i ) Output parameter s : 

i 1ST AT , NPROF, NPLOT : number of axial stations 
between print-out of (i ) station variables, (ii) 

: profile variables and (iii) non-dimensional plots 

XST AT , XRROF, XPLOT : the corresponding axial 
distances (in nozzle diameters) 

The next data cards are read in from LANGL1 and refer to the chemical 
equilibrium constants for the 4 (MR) reactions and the enthalpies for the 
7 (NS) species; these must be in the same order as the auxiliary variables 
are data-typed in Chapter 2 of HAIM. The data are read in at intervals 
of 100°K for temperatures up to 600Q°K, so that for example, HT (2.35) 
refers to the enthalpy of species 2 (oxygen) at 3500°K. 

Further data cards provide the initial profiles in Chapters 4 and 5 
of MAIN. .For those test cases involving step profiles, only the cross- 


20 



stream distances are required together with the values of velocity, temper- 
ature and species concentrations in each free stream. Appropriate profiles 
for the thin mixing layer under consideration are then generated internally.. 
For cases with continuous profiles, the cross stream distances, velocities, 
temperatures and the concentrations (of those species whose presence has 
been previously indicated) are read in and converted to full dimensional 
form. 


The remaining data cards contain alphanumeric data which are used to 
provide headings for the print-out. These are read in from LANGL1 and 
supply information on 


(i) the test case number and description (2 cards) 

(ii) the dependent variables of the calculation 
(10 cards) 

(iii) the chemical rea-ctions assumed (NR cards) 

(iv) the species present (NS cards) and 

(v) their chemical symbols (1 card). 

4 . Discussio n of Sa mp le Predictions 

To examine the general capabilities of CHARilAL, fifteen test cases, 
prescribed by NASA Langley, involving the mixing and combustion of a 
hydrogen jet with various coaxial gas streams have bec-ri computed. The 
computer outputs of these runs have been forwarded separately to the con- 
tracting agency. In this section we examine various features of the num- 
erical predictions, the main attention being given to test cases 1 and 4, 
for which the prescribed initial profiles of velocity and temperature are 
given in Table 2. 

Jet velocity and temperature profiles for Case 1 at four downstream 
stations are shown in Figures 3 and 4. The wake region of the profile 
arises from the wall boundary layers that are present on wall of the 
hydrogen pipe. As the shear f 1 ov; develops downstream tiie jet spreads and 
the velocity level falls. Notice that the jet region is still present at 
x/D = 50*. From the temperature profiles shown in Figure 4 the region of 


In this case ttierc is a slight momentum deficit in the shear flow so, 
far enough downstream, the velocity within the shear flow would be 
everywhere less than the free stream velocity. 
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combustion is clearly evident through the sharp peak 'in the temperature 
level. It is seen that the maximum temperature retrains virtually constant 
as the shear flew develops downstream. 

The decay of mass fraction m ;is shown in Figure 5. Comparison is 

H 2 

drawn between predictions for Case 1 and Case 3 in which initial conditions 
are similar to Case 1 except that the surrounding gas is nitrogen (which 
is non-reacting) rather than air. Also shown is a prediction from reference 
[2] obtained, with the same model of turbulence as that used here, for a 
non-reacting hydrogen air jet with velocity ratio similar to Case 1 but 
where the total temperature of the external stream was only 300°K; the 
density of the air is thus approximately 3 times that for Case 1. (The 
prediction from [2] was in satisfactory agreement with the experimental 
data of Eggers [7], which for clarity are omitted from Figure 5). Clearly, 
because of the smaller density of the external stream, the rate of dilution 
of the hydrogen jet is appreciably slower for the present test cases than 
in the reference [2] computation. The decay rate is seen to be faster for 
the hydrogen/nitrogen mixing than for the hydrogen/air jet. This behavior 
is attributable to the fact that combustion .does not take place in Case 3, 
hence temperatures within the jet are lower and density correspondingly 
higher (if reaction is suppressed in test Case,! the rate of decay of m^ 
is nearly identical with that of Case 3). 2 

Distributions of some of the important dependent variables affecting 
the flow development are shown in Figures 6-9. The station selected is 20 
jet diameters downstream from the exit (through the profile shapes vary only 
slightly with streamwise position). It is seen from Figure 6 that the tur- 
bulence energy k reaches its maximum value at the axis. This is generally 
a feature of axi symmetric jets and is in contrast with the behavior of 
plane jets where the peak energy level occurs near the position of maxi- 
mum shear stress. (The different behavior is attributable to the more’ 
rapid axial decay that occurs in the round jet). Profiles of mean and mean- 
sguare-fluctuation levels of hydrogen element are shown in Figure 7. It is , 
noted that the rapid decrease in the level of fluctuations in the outer 
part of the. jet is the reason that the f + line in Figure 10 lies much 

closer than does f to the stoichiometric line. 

The calculated distributions of radicals and atomic species are shown 
in Figure 9. Their level is very sensitive to temperature and this is why 
a logarithmic scale is adopted for the ordinate, the peak- calculated con- 
centration >f OH is approximately 1% by mass, about ten times as large as 
the maximum' mass fractions of 0 and H. The corresponding distributions of 
molecular concentrations (H^, and H^O) are shown in Figure 10. The 

figure includes predictions for "case 8" as well as Case 1. Initial condi- 
tions for the former differ from Case-1 principally through the presence of 
appreciable water vapor in the external stream. It is seen that, as a 
result, the level of HgO is higher throughout the jet than in Case 1 and 
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the level of N 2 and Op correspondingly lower. The hydrogen profiles for 
the two cases are virtually identical. 

Turning now to the confined flows, it must be said at the outset 
that, for the conditions prescribed in the test cases, the rate of spread 
of the jet was so slow that in no case do the predictions show the jet 
' having reached the pipe wall by the downstream end of the field of compu- 
tation. There are, nevertheless some effects of the flow confinement 
though these are mainly seen in the velocity development (Figure 5 shows 
that the centreline mass fraction of H 2 is only slightly diminished by the 

presence of the pipe v/a 11). Because the flow is supersonic, the exothermic 
reaction leads to a rise in static pressure with distance along the duco. 
This adverse pressure gradient causes a re duction in the level of velocity 
as the jet develops along the duct; velocity profiles at three stations are 
shown in Figure 3. The variation of centre line velocity with x is shown 
in Figure 11 for four test cases. Note that, for Case 6 in which confined 
streams of hydrogen and nitrogen mix without chemical reaction, the varia- 
tion of'U is virtually identical with the corresponding unconfined flow. 
Case 3. **■ 

To convey an impression of the capabilities of CHARNAL for predicting 
both free and wall flows. Figure 12 shows predictions obtained for the mix- 
ing of subsonic f^/air streams. In order to achieve a faster rate of 
spread than in any of the prescribed test cases the initial core jet vel- 
ocity is set to five times that of the surrounding stream; moreover the 
diameter of the confining pipe has been reduced to 49mm. Figure 12 shows 
the development of the velocity and temperature profiles along the duct. 

The jet reaches the pipe wall about 33 jet diameters from the start and 
by 40 diameters downstream the velocity profile looks generally like that 
found in pipe flow. By 75 diameters combustion is virtually complete arid 
the temperature is uniform except in the vicinity of the wall. 

5 . Suggestions for Further Extensions and Refinements of CHARNAL 

CHARNAL provides numerical predictions of how combustion will proceed 
in reacting hydrogen jets. It ought to be emphasized, however, that this 
behovic-r will not necessarily coincide - nor, sometimes, even approximate 
to - the actual development of the flow. Validation and comparison of 
the predicted behavior with available experimental data ought therefore 
to precede any use of the program for detailed design calculations. With- 
out in any way preempting the outcome of such a series of comparisons, it 
may be helpful to indicate briefly some areas where improvements to the 
present version could be made. These are outlined briefly below. 

( i ) Combustion model 


While the equilibrium combustion model embodied in CHARNAL may 
be adequate in high temperature flows, it evidently leads to unrealistically 
fast rates of combustion at low temperatures. The argument in favor of 
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assuming chemical equilibrium is that the chemical constituents may then be 
determined from solution of algebraic rather than/differential equations 
with corresponding savings in computer time. The savings however are rela- 
tively modest because the algebraic set of equations are significantly non- 
linear and require iteration at each node to solve. 

A further important point is that only with a differential reaction 
scheme can one take proper account of the role of /Species fluctuations on 
the progress of the chemical reaction. It is thus recommended that in the 
future CHARNAL be extended to include a finite-rate reaction model. At the 
same time consideration should be given to the inclusion of a more elabo- 
rate reaction model than that implied by equations (2.3-1 - 2.3-4). i 

Turning now to the turbulence-model, it must be said that there are a 
number of areas of uncertainty. The Proceedings of the Langley Conference 
on Free/ Shear Flows (of which reference [2] forms a small part) suggests 
that ttie Mach number of the fluid may exert an influence on the shear flow 
develbpment-at any rate for flows, such as the mixing layer, where turbulence 
levels are high. It is also known (see Reference [2]), tnat the quantity , 

C , although assumed constant in the present work, increases, appreciably 

u ,-l 

when the average level of turbulence energy production at any station be- 
comes small compared with the dissipation rate. Mow it happens that in 
many of the test cases examined by CHARNAL the momentum excess due to the 
hydrogen jet moving faster than the external stream i s just about balanced 
by the momentum deficient of the boundary layers on the nozzle walls. 

This situation leads to a rapid decrease in rate of energy production with 
x. 


If these circumstances are the ones that prevail in the majority of 
tests to which CHARNAL is to be put it would therefore be desirable to 
replace the present constant c p by the elaborate dependency on energy 
production: dissipation rates presented in [2] (associated with the model 

designated ke2 in that reference). Such a change would produce a somewhat 
faster rate of spread of the jet than does- the present model from about 
15 diameters downstream of the jet exit to the point at which the jet 
reaches the pipe wall (if present). There is in addition the possibility 
that combustion will affect in some way the turbulence transport. There 
seems no conclusive evidence on this point yet but this perhaps mainly 
reflects the lack of sufficiently well documented experimental data. 

A separate, albeit related, point to those discussed in the above 
paragraphs is the importance of accurate upstream profiles. An abiding 
feature of weak shear flows such as those tested by CHARNAL is their in- 
ability to 'forget' the nature of their origin. For example the far field 
behavior of axisymmetric wakes depends crucially on the shape of the wake- 
generating object. It follows, therefore, that in order to obtain reliable 
predictions of the present hydrogen/air jets it is necessary that the up- 
stream values of the mean velocity and turbulence quantities be known 
accurately. However, no information was provided (n the profiles cf 
turbulence quantities and so "best estimates" had to be made estimates. 
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however, that may give the wrong levels of k or e by a factor of two or 
three. It would certainly be helpful in clarifying the degree of realism 
provided by the present model of turbulence if a few more experiments could 
be performed for conditions .similar to those examined by Eggers [7] in 
which especial attention was given to a full documentation of the upstream 
flow conditions. 
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7, Nomenclature 


Cjj ’ C l» ^2» Cg]» Cg2 

c . 

PJ 


f 

g 

h 


constant coefficients appearing ;in turbulence 
model . 

mean specific heat of species j. 

diameter of jet nozzle. 

mass fraction of elemental hydrogen, 
mean square fluctuations in f. 
enthalpy. 
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stagnation enthalpy. 

difference between enthalpy. of species at 
298. 15°k and 0°K. 

heat of formation of species j. 

kinetic energy of turbulence. 

3 / 2 / 

dissipation length scale k jc 
mixing length, 
mass fraction of species j. 
static pressure. 

radius (measured from axis of jet), 
radius of pipe enclosing jets, 
temperature. 

local streamwise velocity. 

local turbulent shear stress. 

velocity normal to axis of jet. 

molecular weight of species j. 
distance along center axis 
mass fraction of elemental oxygen. 

effective widths of shear flow, see Figure 2. 

effective turbulent flow transport coefficient 
(subscript denotes diffused quantity). 

rate of turbulence energy dissipation. 

turbulent viscosity. 

density of mixture. 



0 

V 

* 

w 

Subscripts 

t ; 

D 

E 

i 

I 

ref 

U 


effective Prandtl /Schmidt (subscript denotes 
diffused quantity). 

any of the primary dependent variables. 

stream function defined by equations (3.1-2). 

dimensionless stream function defined by equation 
(3.1-1). 

conditions on the jet axis, 
downstream. 

external edge of shear flow, 
value of quantity at initial station, 
internal edge of shear flow, 
reference value of quantity, 
upstream. 

conditions prevailing beyond the outer edge 
of shear flow. 
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Flow Chart of Ck/uuMAL Program. 
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Figure 1. The Flow Geometry Treated by CHARNAL 
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At A, U=0. 9 U T + 0.1 U 
At B, U=U • 

At C, U=0.9 U £ + 0.1 U 

y^ = Characteristic Shear Width for Inner Region of Flow (Axis + B) 
y W2 * Characteristic Shear Width for Outer Region of Flow X (B ->Free Stream) 



a. Continuous Velocity Profile 



Idealized Step 
Profile 



At A, U=0. 9 Uj + 0.1 U E 


At B, U=0.9 U E + 0.1 Uj 

y w ^ = Characteristic Shear Width for Thin Mixing Layer 

b. Step-Change Velocity Profile 


Figure 2 . Initial Velocity Profiles and Characteristic 

Flow Widths 
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Figure 3. Velocity Profiles 
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Velocity and Temperature Profiles in Ducted Subsonic Jet. Initial Profiles 
(a) x/D i = 0, (b) x/D - = 20, (c) x/D, = 40, (d) x/D, = 75 
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Appendix 1. Listing o f F ortran Variables 

Listed below are seme of the most important Fortran variables used 
in the computer program. To prevent an excessive list, attention has 
been confined to those variables which are neither (a) self-evident in 
meaning, (b) indicated by comment cards. in the program, (c) discussed 
in Section 3.5 or (d) connected solely with the finite difference' "formu- 
lation of the equations. Arrays are shown with the maximum numerical 
values of their subscripts. 

Fortran Symbo l Significance 

: constant of proportionality in near-wall mixing length 

formulation, f m = icy^. 

: constant of proportionality in free shear layer mixing 

length formulation, = Ay q. 

L 

: Mach number, u/(yRT) • 

: constants in the dissipation rate equation 


: conversion factor from calories/gram, mole to 

joules/kg. 

’ pressure x molecular weight of mixture 

: constants in mean sqaare concentration fluctuations 

equation. 

: local centreline mass fraction of molecular 

hydrogen * initial value. 

: constant in Prandtl -Kolmogorov viscosity formulation 

\ UJ = Cppk 2 / e 

: constant of proportionality, (= C^) 

CPBAR (7,60) : mean specific heats of chemical species at intervals 

of 100°K. 

CPMN (10) : mean specific heats of chemical species at upstream 

temperatures. 

CSALFA : cosine of angle of inclination of streamlines to 

axis of symmetry. 


AK 

ALMG 

AMACH 

s 

CE1 

CE2 „ 
CFA 

CFAC 

CGI 

CG2 

* 

CLINE 

CMU 

CONST 
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Fortran Symbol 

Significance 

RPDX : 

axial pressure gradient, 3p/3x 

DUDY (40) . : 

mean velocity gradient, 3u/3y 

OX : 

forward step size 


DXLII'i 

; 

limit of entrainment rates (bound to size of forward 
step). 

DYHA 

"• 

growth parameter of jet, dy. ^/dx (see also YHA) 

DDYMA 

: 

k (<, "0.5 /dx > 

EliUL (40) 

• 

laminar viscosity, y 

EilUT (40) 

• 

• 

turbulent viscosity, y-j- 

ENTH (40) 

• 

• 

static enthalpy, h 

F (10,40) 

• 

• 

dependent variables (see Chapter 2 of MAIN) 

FA1 

• 

constant for determination of free-stream 
turbulent kinetic energy and profile of mean 
square concentration fluctuations 

FA2 

; 

constant of proportionality in a free shear layer 
dissipation length scale formulation, t e .=xyQ 

PACI 

j 

location where f-g 1 * = f stojch 

FACE 


location where f*g h = f sto i c h 

FACM 

• 

location where f = f stoic|l 

FRA 

: 

constant of proportionality between forward step size 
and width of shear layer 

FS (10,40) 


auxiliary variables (see Chapter 2 of MAIN) 

FSTOICH 


stoichiometric fuel composition 

GAM (40) 


effective diffusion coefficient 
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Fortran Symbol 
GAMMA 
GASCON 
HO (10) 

I AX 
IEND 

INDE (10) 

INDI (10) 

IPD 

IPROF 

I STAR (40) 

ISTEP 

J1 

J2 

KRAD 

N 

Nfil 
NP1 
OFAC 
OM (40) 

P (40) 

PEI 


s. Significance 

: ratio of principal specific heats, Cp/Cv 

: universal gas constant 

: the, quantities 4h fj2g8 j 5 -, K -(h 298 15 „ K - h Q0l( ) 

(see Section 3.3-4) 

: axial step number corresponding to end of potential 

core region 

: axial step corresponding to, point where jet edge reaches 

duct wall 

: indices denoting nature of boundary conditions 

= 1 (<j> stated), = 2 x (gradient <p stated) 

: turbulence model parameter; = 3 (use of plane flow 

constants), = 1 or 2 [see Launder et al (2)] 

: profile index; = 1 (normalized values), = 2 (full 

dimensional values) 

: number of iteration cycles allowed for convergence at 

each node 

: axial step number 

: 'pccies with largest concentration in equation (3.3-16) 

: species with largest concentration in equation (3.3-17) 

: flow geometry index; = 1 (plane), = 2 (asixymmetric) 

: number of grid nodes, n 

: n-1 

: n+1 

: ratio of oxygen element/nitrogen in outer stream 

: non-dimensional stream function 

: static pressure 

: change in stream function across flow, ^ 
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Fortran Symbols 
PR (10) 

PRDRP 

PRT (10) 

PS.IE 

PSII 

PSIR 

QE 

R (40) 

RC (5,60) 


RCON (10) 
REXD ■ 

RFLOM 
RJE (10) 
RJI (10) 

RJTE (10) 
ROT I (10) 
RME 
RMI 


Significance 

molecular Prandtl /Schmidt numbers 
pressure drop, - p 

turbulent Prandtl /Schmidt numbers 

value of stream function at external boundary ^ 

value of stream function at internal boundary vpj 

value of stream function at duct wall 

wall heat transfer rate 

distance from axis of symmetry ( = Tj + Y cos «) 

chemical equilibrium constants for each reaction at 
intervals of 100°l< 

[N.B. RC (5,1) contains constants for the global 
reaction, H 2 + %0 2 = H 2 O] 

chemical equilibrium constants at upstream conditions 

j 

dimensionless excess radius in pressure gradient 
formulation for confined flows 

flow radius (see REXD) 

diffusive flux at external boundary 

diffusive flux at internal boundary 

total flux (convective and diffusive) at external boundary 

total flux (convective & diffusive) at internal boundary 

mass entrainment rate (radius mass flow rate) 

mass entrainment rate (radius mass flow rate) 
at internal boundary 


ROUBAR 


: mass flux for confined flows ( = mass flow rate/ 

X-sectional area of duct) 


RTW (40) 


turbulence frequency (dissipation rate/ turbulent 
kinetic energy), e/k 
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Translation 


Fortran Symbols 
SCH (40) 

SD (40) 

SU (40) 

TAUE 

TBAR 

TUNE 


Prandtl /Schmidt number array for near-wall flows 
storage location for "downstream" source terms 
storage location for "upstream" source terms 
shear stress at wall 

bulk temperature ( / r UTdy / / r Udy) 

duct duct 

local centreline temperature/initial value 


UBAR 

UEIM 


: bulk velocity ( = / r U 2 dy / / r Udy) 

duct duct 

: initial centreline velocity excess (u.. - u ). 

X <*> 1 


ULIi'.E 


ULIMI 


ULINE 


limit on entrainment at external boundary (to prevent 
formation of profiles with 'long tails') 

limit on entrainment at external boundary 
(to prevent formation of profiles with 'long tails') 

local centreline velocity excess (u^ - u co )/(u^ - uj^ 


vrnx 

Hi i (10) 
X 

XAX 

XD 

XU 

XUD 

Y (40) 
YHA 


: reciprocal of molecular weight 

: molecular weight 

: concentration of oxygen element 

: length of potential core (in nozzle diameters) 

: axial distance at downstream station i j 

: axial distance at upstream station 

: upstream axial distance/ jet nozzle diameter 

: transverse distance from internal flow boundary 

: jet half-width, y 0 
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PROGRAM CHAPN«L (INPUT* OUTPUT, TAPE5=INPUT. 7 APEp=OUTPUT ) 

HI HI <1 Hi HI \l HI Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi 

C 

c spalding and patpmkak passa programme for pound ary layer flows 
C MOD 1 F I EO BY a. morse for PREDICTION OF CONE INED/r *<EE H YOROGEN-4 I - 

c mixtures fo« langley RfSfiarch center, u.s.a. june r-r/3 

c 

c HHMHUi Hi Hi <» Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi ii Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi 

DIMENSION THETA (AO) 

COMfloW/C>ENRAL/ACC. CSALFA ,OPDX .DX.ENTh (40 ) ,F ( 10 ,40) , FS (10,40) , GA'i j, 
1 bAMU.GAM ( a o ) • I , IF. IN, I NOR (10) , I NO 1(10), I ST r.P , I UTR AP ,T TEST , KEX , K I ■ . , 
?KPAD,M»NE0,NP1 ,NM1 ,0" (40) , OHIO (40) ,BOm (40) , K AST ( 1 0 ) , PE I , PS I E * PS 1. 1 , 

3 P ( 4 0 ) , K ( 4 o ) » R h 0 ( 4 0 ) , R J E ( 1 i> ) , R J I (10) *PMF. , WM I , SD ( 4 0 ) , SU ( 40 ) , WM ( 1 0 > , 
4XD,X|i, Y (40 ) t Yl) IF (40) , YE, Yl ,RJTE (10) , RJT 1(10) 

C 0 M M o N / C R a D / R F L 0 W , R E X D , K A S E 

COMI-iOM/Cjs/JtJ, JK, JD, Jhs* JA, JG, JTF , JLE, JUV, JB 
COMMnM/CFF/JHP, J02 , JOH , JH20 , JH, JO , JN2 

OOMMOM/CTpHH/AK , ALMG • CMIJ , CMOIN, CE1 , CE2 , CG 1 , CG2 . CP I T « IHJi)Y (40 ) , 
IDUDYSO (40) ,KMUL (40) ,r NUT (40) , IPD,PR(10) ,PRT ( lw ) ♦PTW(40) 

EOMMpM/OUTP T / I a X , I END , Y I N , YOHT , YihA , UYH A , DUYhA , RDUCT , KP , KS , K-T , 

1 KASFT'0*XSTAT ( 2 0 ) , XPROF (PO) ,XPLCT (PO) , PRESS 
COMMpM/CwwF/YR (p) , J RuF <2) , EWALL * H 

■COMMoM/TPLOTH/XTAXiS.XtPI.OT (40) ♦ YT.AXES (10) , YTPLOT (10,40) , 

) YTMAx ( 10) , YTSYMP (10) , OUT (40) , IPROF 
rOMHOM/CPROP/lO? , IISI 2, I H?0 , OF AC , NR . NS , WOO , wRO , WSO , WTO , G A SCON , G AMP A 

CHiHi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi Hi HI HI HIHI Hi HI HI HHMI HI Hi HIHlHi Hi Hi HI HI HI HI Hi 

CHAPTER 1 

C4HHHHiP/,RAt'FTER<; AND CONTROL INDICES 

R E A 0 f S » 1 1 5 0 ) KASRNO ,N,NEO *KASE ,NSTAT ,NPROF,NPLOT , 1PD* K IN, KEX , KOi'JF I 
IN, IOp, IMP', 1020, NR, NS, IljNF 
Rf : AD (S» J 1*52) Xu, XJJLAST » YIN,YOUT»M[)UCT*PReSStTAtUlW 
IF { ItiMF.EO.O) RpAD ( 5 , 5 0 1 ) Y W 1 , U R 1 , 1 W 1 
IF (K/.sE.E^.^) RF AO ( 5 , 500 ) TDUCT 
PPESSsPRFfS^I ,E*05 
X( i = X| ]H> Y OUT 
X(ILAST = XULASTHIY0UT 
YJN = YTN»Y,nUT 
DPDX = 0. 

READ f 1 ? *1152) (XsTAT (I ) »I = 1 *MSTAT) 

RFAD > S * 1 ISP) (XPROF ( I) , 1 = 1 ,NPROE ) 

W L A D ( S * I 1 S P I (XPLOT (I) ,1 = 1* NPLOT ) 

XRROF(NPRpF+1)=XSTAT (NSTAT+1)cXPlOT (NPLOT+1)=1«E+30 / 

CALL S T R J D 0 

I ax=j nuooo 

TFNDsl 000 00 

IF (KFX ,EO. 1 ) I E N D = 0 

LASTFP=20ft0 

I PROF a2 

K R = .1 

KS=1 

KT = ) 

ML IMls .01 
ui- 1 fir= -oi 
DXLlf's.O'j 
FR A = , OR 

AF A=.?S „ „ 

c «**largfr forward step for test case no,i6 to cover longer distance 
JF (KasENO.EO. 16) FHA=.125 
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On 1 <\ J = 1 . N £ 0 
.14 I f-J 0 E ( J) =1 
I f'i D E '( JA) = 2 
imjF \ \) =IR!jF (?) =0 
1VC~=2 

C<K><> <MH> ##«#«#»## »#»<(««»«#»«# tt ««»»« »*#<(■#«■ 
C H A P T f R 2 

chhmhm»sele;ctxon of dependent variables 

DATA JU, Jk, JU. JHS» JA, Jq, JTE» JLE* JU V . JB/ 1 , 2 , 3 , 4 t 5 , 6 , 7 , 8 , 9 , 10/ 

C 

C F ( JU, I ) . , . . AXI AL VELOCITY 

C F <JK, j >. ...KII-ietic- ENERGY OF TURBULENCE 
C F(Jf),i) ....dissipation RATE 

C F (JUS. I ).. .STAGNATION ENTHALPY 

C F ( JA, T) ... .HASS FRACTION of HYCRoGEN SPECIES (H2,H,h20,OH) 

C F(JG,I) ...MEAN SQUARE CONCENTRATION FLUCTUATIONS (OF H ELEMENT) 

C F(jTr< I ) ... absolute TEMRFRATURE 

C F ( J L F • I ) ... DISSIPATION LENGTH SCALE 

C F(JUV»I) ...REYNOLDS STRESS CORRELATION (UV) 

C F(JB,I) ....MASS FRACTION OF SPECIES B 

C 

c ^auxiliary variables 

c 

C FS (JHP* I ). .CONCENTRATION OF H? 

C FS (JH, I ).. .CONCENTRATION Of II 

C FS (JOH* I )., CONCENTRATION OF OH 

c fs (J ppo, i j .concentration of H 20 

C FS (JOP* I ). .CONCENTRATION OF 02 

C FS (JOtl ),. .CONCENTRATION OF 0 

C FS (JN?» I ). .CONCENTRATION OF N2 

C 

DATA JH?, JO<S JOH , JH?Ot JH, JO, JN2/1 ,2, 3, A ,5,6,7/ 

C <HHH!OLECULAR WEIGHTS (AND THE I R RATIOS) 

DATA ( wm { j) * J=l, 7) /2.0 16, 32. 0» 17.008, 18.016,1.008, 16.0*28.016/ 
WKo=WM ( JO) /WM ( JOH) 
w T 0 = 1 • *• W ,7 n 
W S 0 = W 1 1 ( J 0 ) / W M ( J H 2 0 ) 

WA- 0=1 ,-WSO 

C-IHt <HH> IHHHUHttHHUt <> QiUHHHt «#«««#« 

CHAPTFR 3 y 

C<HHMHtcONSTA<J r s 

c HHHH 1 ATEPTAL oONSTANTS (S.I. UNITS) 

GASC0N=B314. 

GAMMAsI.4 

C «##C'HF |v ' I CAL FOULIRRIUM CONSTANTS AND ENTHALPIF.S 
CALL LANGL1 

C Tmh>TI.IRBi'LENCf CONSTANTS 

A K = . 4 0 

ALMGs.l 1 

CMUIN=«09 

CONST=l ./SORT (CMUIN) 
rfijafMUiN 

cri^i .43 J 

r.F?=i .v 2 
001 = 2.0 
C02 = ( ". o 

C SPECIFICATION OF PRANDTL /SCUM I OT NUMBERS 

PPT<JU>=1.') 


50 



PPT (JK> =1 . 0 
PUT (JO) =1 . 3' 

PPT ( JHS ) a , 7 
PRT(JA>=.7 
PRT(JR>=.7 
PRT(JG)=.7 
DO 4nJ J=1,N£Q 
401 PR ( J) el • 0_ 

FA l =4 . £-04 
FA2=,R75 

fviall= 9 » 

H=. 9 
CPITeO. 

C 44444 44444 4# {>«4H»44'i)' 44444 444 44 44444444444444 4444444 44444444444444 
CHAPTER 4 

C*<hhh»sPECiFICaTION OF GEOMETRY 
KPAUsP 
CSALF4=1 . 

R(l)cYlM 

R E A 0 ( 5 » 5 0 0 ) (Y (I) *1 = 1* NP 1 ) 

00 40 1 = 1 .NH1 
iF(IUNP.Fo.l) Y<n=.l«Y(I) 

40 Y(I)sY<I)^YOiJT-YIN 5 

IF (KFAU.Eo.2) GO TO 163 
DO ip? I el » NP 1 
Ifl2 Hills!. 

GO TO 161 

1R3 DO 204 1=1, NP l 

POO R< I >=P<1) *Y < I ) #CS ALFA 

C44444444444444 44 4 4 4444 444444444444444444444444444444444444444444 
CHAPTER 5 

0*4*4# INITIAL CONDITIONS 
101 DO 4S0 1 = 1, NP1 
450 P ( I ) ePRESS 

IF(ImWF.EO.I) GO TO 723 
C 4440.ONTTNOOl.ls PROFILES 
YW1=YW1 # Y()UT 

YWl=ywl-YjN 
l.|W 1 e(l W 1 4|J j N 
IW2bjw1*i 

READ { 5 > So 0 ) (F(JUfl ) ,I = 1,NP1) 

R F A D 1 5 » 5 0 ii ) (F ( JTE , I ) » 1 = 1 »NP1 > 

RFAD (5 » So 0 ) {Fs<JH?,D .I = 1,NP1) 

IF(Io?.fc.Q.l) PEAO(5,50 o) ( F S ( JOE . I ) , 1 = 1 *MP 1 ) 

IF(INP.EQ.I) «EAO(5,SOO) (FS ( JN2, 1 ) , 1 = 1 ,NPl ) 

TF(lH?D.E0.1) READ(5,500) (FS ( JH20, I ) , 1=1 .NP1 ) 

C 444MASS FRACTION OF HYDROGEN ELEMENT 
DO InS 1=1, NP1 

105 F (JA'. I )=Fs (DHP, I) + FS (JH, I) ♦WQO#FS ( JH20, I ) ♦ WT0«FS ( JOH , I ) 

C 444RATIO DF oXYGFn ELEMENT TO NITROGEN IN OUTER STREAM 
OFACsFS (Jr2»NP] ) + WS0«F S ( JM20 , NP 1 ) 

OFACioPAC/ (OFAc + FS ( JN2.NP1 ) ) 

C 444C0NVFR3I0N OF INITIAL PROFILES 
DO 4 s I = r.NPl 
F ( JU, T. > =F '( JO, J ) 4U IN 
45 F(JTF.1)=F(DTE,I)4TA 

C ##4TURm.H.tNT KINETIC ENERGIES BY MIXING LENGTH HYPOTHESIS 
On lofl« T=?,N 
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Dl-'DY f I ) = (F ( JU, 1*1 ) -F { JU, 1-1 H /TYTl*l )-Y <1-1 l ) 
in An IIUOY { T ) =I) 110 Y ( 1 ) *miUY ( I) 

Ol'OY ( ) ) =n|i[)Y (DPI ) =f). 

c <nn»Df-FIr\'F CHARACTERISTIC shear layer widths 
MR£Fb//«F"(JH.1 > . 1 *U‘J( 1 
Du lp62 1 = 1 » NP 1 

IF (F (ju, n .LT.URfcE) GO TO 1063 
1062 C.ONTlMUc 

10 63 yv? = Y fl-l) ♦ (Y(l)-Y(I-l ) ) ( UREF-F ( JU, 1-1 ) J / (F ( JU, I ) -F C JU, 1-1 )) 

YWP=YWl-YW? 

Y'J2 = ALH<,*yW* 

JF (KAStHO.EU. 16) GO TO 8702 
iiREFs.. Hmiw1 + .9»F (JU,NP1 ) . 

DO lp'>4 I = I W i * N P 1 
IF (F (,JU. l) ,GT .UREf) GO TO 1065 
1064 CON TT HUE 

10 65 YW3 = Y (!-!)♦ C Y ( I ) -Y ( I- 1 > ) <> (UREF-F { JU* 1-1 ) ) / (F < JU, J)-F < JU* 1-1 ) ) 
Y w 3 = y U 3 - y i.i 1 
Yw3 = /! |J'10 «yW3 

8702 IF (KaSEMO.EU. ) 6) Yw3=YW2 
DO 10 66 I = 1,1 ui 

1066 F(JP,T)aCcNST#YW2«YW2*0UDY{I) 

DO 1067 I = I W 2 ♦ N P 1 

1067 F(JK,l)=C0M8T«YVJ3*YW3*r)UDY(I) 

Do l 068 1 = 1 , DPI 
•akmim=kai»f < ju, r >*f (JU,d 

1068 F ( JK, T ) =AMAAl (F < JK, I ) ,AKMlN) 

c dissipation hates 

Y'«‘?=.Y'.J2/aLM.G*Fa2 ...... . ' 

Yw3=ywJ/ALMG»F A2 
do 1 n6‘> i = l, l ; <i 

1069 F(JD,T)=FfJK,.r>#SORT(F(JK»I))/YWP 
DO 1 0 7 0 I = IW2.NP1 

1070 F(JD,I)sF(JK,I)*SGRT(F(JK,I) )/YW3 
00 To A 17 

723 roNTjMUE 
C <hm»sTEP profiles 

«FAi.)<5*5on> UI,uE,TI.TF,02£,AN2E,H20e 

HP 1 Dp = /'JP i /2 

c ihh > con struct linear velocity profile 

DO 11 Ial’.NPl 
YHAT = Y < I ) /Y (DPI ) 

F ( JO, T ) =iJT*YHAT» (UE-Ul) 

IF ( J. .GT.MP1D2) GO TO 454 
F ( JTf , I ) sTI 
F S ( JF? , I ) = 1 • 

Fs ( JpP » I ) =0 . 

FS { Jf 12 » I ) =0 . 

FS { JH20, I ) =0. 

GO To 11 
454 F { Jig- ; I ) =TE 
FS (JH2* I ) =0. 

FS ( J02, I ) = 0 F F 
F S { JN2 » l ) =AH?E 
FS ( JHZO, I ) =HEOE 

11 F< JA, I ) =FS ( JH2,I) + F S ( J H , I ) ♦WOO«FS I JHPO, I ) ♦ W TO«FS ( JOH , I ) 
OFACcFS { jp.2 , NP 1 ) +WS0»FS( JH20,iMPl) 

OFACsOFAC/ (0FAc*FS( JMZ,KPU ) 
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alnth=»p«y (NP1) 

IIGHA0=<UE-UI>/Y(I'P1) 

LIPRAOcUGR aD^ALWG^ALNTH 
UGRAn=COM S T«UGRADouGRAD 
C **#TURtUiLkNT KINETIC ENERGIES 
DO 472 1=2. N 
A72 F (JK, l)=upRAI) 

DO 473 T=I ,NMl 
AKMlH=EAlttF(JU,n«F(JU,I) 

473 IF (F ( JK, 1 ) .UT.AKMIN) F (JK.I)=AKMIN ...... . 

ALNTH=FA2»ALNTH f _ 

C ##«*[> I SS I PATinN RATES 

DO 474 1 = 1 ,NP1 . . / ’ 

474 F ( JD t I )=F > jK, I )tt$ORT (F ( JK,I ) j/ALNTH 

C #<n>«;TAG|UA'riON! ENTHALPIES 

417 DO 4c 1 = 1 NP1 • ; ’ . 

CALL LAMc,L2 • . • . 

49 F( ! JH S ,I)=fNTH(I) *,S»F<JU,1>»“2 + F(JK.I) 

C «*«MEAN SQUARE CONCENTRATION FLUCTUATIONS (OF H ELEMENT) 

DO 12 1=1, NP1 

12 F('JG',I)sFAl*F(JA»I)*F(JA,I) 

Q iHHH> » Q <*■<*# ■"** ■<*' f # ### 1HHH* p <> V * 
CHAPTFR 6 

C<HH.4»nMEGA DISTRIBUTION 
DO So 1=1, NP1 

VMIXcO. 


DO 107 J=1«NS 

107 VMIX = VMIX + FS ( J. I ) / W M ( J ) 

RHO ( j ) aP ( t ) /GASC.ON/VtllX/F (UTE. I ) 

50 THETa(I)=RHO(I)»F (JU r Il*R(I> 

DO 5p5 I=2,NP1 

ZFTAs.b* (THE FA ( I) ♦ THETAU-1) ) * (.Y (I) - Y (1-1 ) > 
505 ON ( I ) =0M ( j- 1 ) + ZFT A 

PSII=RHO(l)«F(JU,l) tt YlN 
IF (KrAD.fo.2) PSII=.5«PSII0YIN 

PSlEsOM(NPI) *PSIJ 
PF I=0 M (NP1 ) 

00 5o 4 I =j » NP 1 
504 OMII)sOM(i)/PEl 
CALL S T R 1 0 1 




CHAPTER 7 

C*#tt##THERHODYN/iMlC AND TURBULENCE PROPERTIES (START OF MAIN LOOP) . 
100 CONTINUE 

TF ( IsTEP.cE. IAX) KJN=3 
IF (I-sTEP.pE.IEND) K F. X = 1 
C #*#LIMlT5 ON SPECIES MASS FRACTION 
FAMll.isl »E-30 
F AMAX = 1 • 

DO 455 1=1 »nP1 

FfJA, T)= amIN1 (F(JA,I) »FA«AX) 

455 F(JA.T)=AmAX1 (F(JA,I) ,FAMIN) 

C ##«L IMlTS ON TUKHijLFNCE PROPERTIES 
FGMlftcl .F -30 
FKMlNal.F-30 


DO 600 1=1, UP1 

F ( JG', I ) = amAX 1 (F (JG,I) tFGMlN) 
FKMAx = F ( Jl'» I ) *F ( JU, I ) 
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F ( JK, T )=Ah'INl (F ( JK, I) iFK^AX) 
F(JK,T>=AMAX1 (F(JK,I> tFKMlN) 

FDMlNaF ( Jl', I > «SOPT (F { JK» I ) ) /Y (NPl > 

*00 F(JI>.I>=aPaXI (F(JD,I) tFDMlN) 

dp=dpdx*dx 

iF'dpTtP.CO.O) GO TO 425 
c ***LftCAL MAs<; FRACTIONS 
CALL LANGL3 

C <nn»LOC*L PRESSURES AND TEMPERATURES 1 
DO 60l 1=1* NP 1 
P(I)=D(I)+DP 

601 FNTH(T>=F(JHSfI)~,5*F(JU*I)<n»2-F{JK*I> 
CALL LANGL4 
DO 6ft? 1=1* NP 1 

IF{F(JTE,t) .GT.O.) GO TO 1083 
WRITF f6,2702) F ( JTE * I ) » I « ISTEP 
TFlNal 


1 ftfl3 VHIXsft. 

DO lftP J=1*NS 

1 09 VH I Xs V/M I X *FS ( J » I ) /wM ( J) 

*0? RHO ( I ) =P( T) /GASCON/ VM I X/F < JTEtll 
*25 CONTINUE 


IF(lrTN.FQ.l) GO TO 117 

C«« »#*<>«*<»«»***<>* ****** <t< *** <t ' < * tt ** ************** 


CHAPTER 8 

C*«<hm»TRaNsPORT PROPERTIES* PRESSURE GRADIENT AND FORWARD STEP 
IF ( K j N.NE . 2 ) GO TO 522 
IF (KraIj.Eo.2) GO TO 521 
YIN = ps1 1/ (RHO ( 1 ) «F( JU*D) 

GO To 522 

521 YIN = coRT<A0S(2.<>PSlI/<RH0<l)*F(JU»l) ) ) ) 

R(l)=YlM X 

c 4 **$PHEa01NG rate 

5?2 CONTINUE 

CALL YlNT’, ,5,YHA*JU> 

IF ( ISTEP. I-Q.O) GO TO 239 
DYHAa ( YHA-YHALS) / (XU-PXU) 

DDYHA= (OYHA-DYHAV) /DYHA 

239 YHALS=YHA 
()YHAV=OYHfl 
PXU=XU 

238 CALL STRID2 

IF (ISTEP. NE.O) GO TO 1370 
IF (KagE.Nf.2) GO TO 1370 

c **«stream function at duct wall 

PSIR=pSIe* .5*RHq (NPl ) «F ( JUfNPl ) * (RDUCT«*2-R (NPl ) **2) 
ROUHAR=2.«PSIR/ROUCT#«2 
1370 IF(KEX.EQ.I) K'0NFIN = 1 
GO TO (71,72) * KONF IN 
71 .TF (ISTEP, FQ. IEND) RDUCT = R (NPl ) 

IF (IVCE.Ep.2) GO TO 1338 

C #»»PaSSa VERSION OF PRESSURE GRADIENT FOR CONFINED FLOWS 

C •cM»<»cALClil.ATinN OF PRESSURE ADJUSTMENT 

C — VERSION uirTH GRID FILLING THE DUCT 
IF ( ISTEP. FO.O) GC TO 72 
a(>uCt=- 5 « (HiiucT#unnCT-R ( U «R ( 1 > > 
AFL0Wa*5*jrH(NPl)*W(NPl)-R(l)»R(l) > 

PSlOIFsPEi 
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RnijHAR=PSTDlF/ADUCT 
nRAR=0 . 

C>o 73 1=2'.N 

73 |iBAR=llBAR4..5*H() M (r)#F{ju,I J 
roqap=hourah/ubar 

DllDPs-1./ fROUfiAR* (RJE < JIJ)-RJI ( JU) ) #DX/ADUCT/UBAR ) 

UPROpP=R0MBAW/fiAMMA/PfiESS 

DROU=ROUBaR*(AFLOW/ADUcT-1.) 

DPsDRnO/ (ROBAR*OUDP>UDROOP) 
i lFAC=r)P*Dl)DP/UBAR + 1 • 

c #*« ADJUSTMENT of velocity, pressure and density 
DO 7A I = r,NPl 
F ( JU, I ) =F ( JU* I ) *uFAC 
P ( I ) aP ( I ) *DP 

74 PH0(T)*RHO(I)*a.*GAMMA*DP/P(I)) 

DPDXbOPDX+OP/OX 

C #**RE CALCULATION OF DISTANCES 
CALL STRIDE 

oo to f?. 

c *»*(3l'NMIX VERSION of PRESSURE GRADIENT FOR CONFINED FLOWS 
133s sTORKsTHAR 

T Fi A R r a . 

IIRAR = s. 

00 1340 I=?,N 

TBAR=TOAR*.b#bOM ( I) *F ( JTE, I) 

I 34 A IIBAPsIIUAR* .Sorom ( I } »F < JU, I) 

IF ( IsTLP.ftE.IEND) RFL0W=R(NP1) 

TF ( IsTtP.ftE • IFNO) GO To 0719 

TbARa (PSlTttF (JTE,1 ) ♦PE I#T«AR+ ( PS IR-PSIE) #F (UTE ,NP 1 ) >/PSIR 
1.113 ARs (PSj r«F< JU, 1 ) ♦PEI*URAP+ (PSIR-PSIE)»F<JU»NP1) )/PSIP 
IF (ISTEP.fQ.O) 00 TO 7 2 
[)P = -f:nUHAR«UHAR* ( 1 .-STORE/TBAR) 

RFLOWsSQRT (R (NP1 ) <m»2 + 2.« (PS IR-PSIE) / ( RHO ( NP 1 ) #F ( JU,NP1 ) ) ) 

9719 CONTINUE 

IF ( IsTtP.ftE. I END) DP = -R01J0AR»UBAR# (1. -STORE/TBAR) 

REXUs { RFloW-ROijcT ) /RDUftT 
DP=OP-AFAttROUrtAR«UHAR«REXQ 

IF (ISTtP.ftT.IEND) DP=DP-2.*RjE(JU)/R(NPl)*tt2»DX 
DPsOp/ ( 1 . -ROUBARttUBAR/P ( 1 ) ) 

DPDXsDP/DX 
C 4HM»F0 RWaP 0 STEP 
7? DX = FpA<»Y (NPl ) 

I F ( Il)NF . fo . 1 ) [)XsFRA«R (N^l ) 

IF (IsTtP.LT.bO) 1)X = DX«FLOAT(ISTEP+1)/SO. 

IF ( I.STtP.pE* I END. AND. I STEP* LE. I END + 9)- DX = . 1 «DX*FLOAT ( I STEP- 1 END* 1 ) 
IF(ARS«REXD» .GT..005) DX = OX# . 005/ABS < REXO ). 

DXsAPlN-1 (DX»XllLAST-XU) 

X D = X l) + D X 

C ### # <H> tt « «««# » «####» ****# 

CHAPTFR 9 

CHHHHH> a DJUsTMENT OF BOUNDARY CONDITIONS 
C ROUNDARY VALUES ADJUSTED IN STRIDE (3) 

TF (Kji\I.Eo,2) GO TO 95 ■ 

RMIar. 

TF (KpAU.Eft.2) R() ) =0. 

yin=o. 

psi I= n. 

95 IF (KFX.NE.U GO to 196 
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F(JU,nP1)bO. 
rjf ( ja > =o. 

F ( JK ; mP1 ) =0. 

F ( J!.),NP1 ) =0. 

F ( J(5 j !\|F ] ) 3 0 • 

F(JTF,NPi)=.TOUCT 

r=NPi 

CALL LANfiLF 
. F(JHs.Nf»i) =t‘i'vTH(NPn 

CHAPTFR 10 

C<kmmh»TRaNspORT aNO ENTRAINMENT properties 
IP* CALL AUX(JU) 

IF ( I$TE.P,p(£. IEN01 HMF = 0. 

c *#«knti«a tnhent control 

C «*«RFSTPTCr ATTENTION TO VELOCITY AND TEMPERATURE CHANGES 
IF { K I M • M{- . 2 ) &o TO 9 A 
RMI=r (3) !^aM3/y (3) ... 

RMI=RM1/prt ( JA) 

K A T = APS ( (F (JiJ«2)~F (JUtl) ) / { F ( JU » NP 1 ) -F ( JU * 1 ) +1.E-30) ) 
AMGN=AMS( (F<JTE,E)-F(JTE»D ) / (F (JTEiMPI)-F { JTE » 1 ) ♦ 1 . K-30 J ) 
RAT=AMAX1 (RAT , amGN) 

_ TF (RAT.LT.ULINj) RmI=RMI«RAT/ULIMI 
9a IF(KFX.NE.P) GO TO 97 

RMfra.w (NM1 ) *GAmN/ ( Y (NPl ) -Y (NM1 ) > 

RME=RME/PRT ( JA) 

R A T = A 0 S ( (F (JU,N) -F( Jlj.NPU) / ( F < Ju»NPl ) -F ( JU , 1 ) ♦ I . E-30 ) ) 

Alir,U=AMS{ fF (JTE,N)-F{JTE.NPin/(F(JTE,NPl)-F{JTE,l) *l.E-30) ) 

R A T = A M A X 1 ( R A T ? A H G N ) 

TF (RAT.LT.UHME) RmE=RME*RAT/ULIhE 
97 IF ( ( f; M I - R f' £ ) ) X , L T , P K I * D X L I M ) GO TO 96 

OXal)xi.lM«PEl/ (RMI-RME) 

XDaXll + NX 
9 A CONTINUE 

IF (KASfc.EO. 1 ) GO TQ 960 
IF ( ISTEP.pE. IENO) GO TO OaO 

C «*»aOjUsTNFnT OF FORWARD STEP FOR JET TO . REACH DUCT WALL 
TF (« (MH1 ) ,LT, ,999«RDIJCT) GO TO 960 
ieno=tstep+i 
X rNDsXD/YnUT 
'•'PITf (0,5n8) I END *XEND 
9 A 0 CONTINUE 1 

C #*«AOJl:ISTMENT OF FORWARD STEP TO REACH AXIS OF SYMMETRY 
IF (Kin.Ne, 2) GO TO 195 
IF(PSII.GT.KMI#0X) GO TO 195 
OX=PsiI/RwI 
X I ) = X u + 0 x 
IAX=ISTEP*I 
XAXaXD/YOnT 
WRITE (6*5n7> IAX.XaX 

C H A P T T R 11 
C<MMH»«mjTPUT 

195 TF ( ISTEP.NE.O) go TO 193 
CALL OUTPi 
CALL oUTPp 
CALL 0UTP3 
HQ To 113 
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193 XDD=XH/YOl(T 

TF (XUD.LT. XSTAT (KS> ) GO T'< 

KS=K?*i 
CALL; OUTP? 

321 TF jXijn.LT.XPHOF/KP) ) 60 TO 112 
KP=KP+1 
CALL OUTp.3 

11? IF(ISTEP.NE.IAX.AND.lSTEP.NE.IENO) 60 TO 113 
117 CALL oUTp-i 
CALL 0UTP3 

. IF(IFIN.EO.I) GO To 1 0 0 1 

C *«*FoRMAT STATEMENTS 
sno format (7fi o.3) 

SOI F'ORMflT (HFIO'3'U) 

507 FORMAT (/lHO* tt --MlXiNG layer REGION Ends AT ISTEP=*,I5»3X,*L£NGTri o 
IF POTENTIAL CORE =*,F7.2,* DIAMETERS*/) 

500 FORMAT (/1H0»# — JET REACHES DUCT WALL AT I STEP=* * 1 5 , 3X , *DOWNSTRE AM 

l distance =**>/,?.* diameters**/) 
liso format ( )o'r a) 

115? Format ( 1P7E1 1 .3) 

270? Format (/lHOf** NEGATIVE temperature OF*, lPEll. 3. * calculated at 

INODE*. 13, * AT I STEPs* , 1 5 ) 

C ********* * * * * * * * * ******* ******* ******** ********** * * tt * ***************** tt 
CHAPTER 22 

C****«CONTIMUATTON> TERMINATION 

113 IF ( IsTtP.cE.LASTEP.OR.XU.GE.XULAsT) GO TO 1001 
CALL S T R 1 0 3 
GO To 100 

1001 CONTINUE 
STOP 
END 


compiler space 



Sl'RKn'lTlNF AUX(J) 

CPMHnN/GENRAL/ACC*CSAUFA,OPDX,DX,ENTH f40) ,F (10,40) ,FS( 10,40) r GAM3 . 
1GAMN,GAM ( 40 ) * I , IF'IN* INDE (10) , IfslD I (10) , I STEP , I UTR AP , I T EST » KKX , K-IN , 
?l<RA(),N»Nf;o ,NP 1 , NP l , ()M (40 ) »OMD (40) » ROM ( 4 0 ) » K AST ( 1 0 ) » PE I > PS I £ . PS 1 1 , 
•RP (40 ) ,R (40) t Who (4 0) »RJE ( ) 0) , RJI (10) t RME » RMI , SO (40 ) *SU (40) ,WM ( 1 0) , 
4X0,X(U Y (4n) *YDIF(40) . YE , Y I , R JTE ( 1 0 ) »RJTI (10) 

OOMMON/CJs/JU . JK. JO. JHS. JA, JG» JTE. JLE» JUV. JB 
rOMMoN/CFs/JH2, jo2, JOH/JH2Q, JH, JO* JN2 . 

FOMHON/CT|)R:i/AK,ALMG,CMU.CMUIN,CEl.CE2»CGl ,CG2»CRITfDUOY(40) * 

-• 10UDYSO<40) ,EMUL (40) ,EMUT (40) »IPD,PR (10) ,PRT (10) *RTW(40) - 
COMMoN/A(JXST/TOA (40) »U1P 
COMHnW/Clr-MA/SCH (40) 

C<hhV«h> aX I A|, VELOCITY 

. , IF ( J.NE. Ji'/) GO TO 299 

C <hh>L A M I U 4R VtSCOSITY (APPROXIMATE POWER LAW RELATIONSHIPS) 

C ---ignore atomic AND radical concentrations 
DO 60 I = l',NPi 

SUM = FS < JH?»I > *FS < J02. I ) +FS ( JN2, 1 ) *FS ( JH20. I ) 

EMUL (I>=0. 

C —HYDROGEN ( H 2 ) 

EMljL ( I ) =FMUL ( I ) *FS ( JH2* I ) «a.42E-0b« (F ( JTE. I ) /273. ) «»0.675 

C OXYGEN' (OP) 

FMUL ( I ) =EmUL ( I ) +FS ( J02 . I ) *1 .92E-0S# <F ( JTE. I ) /273. ) »» 0.735 
C NITROGEN ( N 2 ) 

FMUL < I ) -EmuL (I ) +FS(JN2, 1 ) <d.66£-05»(F(JTE*I)/273.) **0.690 
C WATER VaPpUR 

EMULr I)=E mUL(I) +FSIJH 20.I) *1 .71E-05* (F ( JTE. I ) /50 0. ) **1 . 050 
69 FMUL ( T ) = ErljL ( I ) /SUM 

C ^^DETERMINATION OF ACCELERATION PARAMETER 
ACC=ARS (ACC-AHS (ACC) ) 

CALL Y I N T ’( . 99.YRS.JU) 

ACC=. G^YHc^ ACC/A R$(F (JU.l ) -F ( JU.NP1 ) ) 

ACC=ACC**n,2 

IF ( IsTEP.FQ.O) GO TO 74 

C ■*<h»FVaLuaTION of INTEGRATED production/dissipation ratio 
IF( IPD.NE. 2) GO TO 73 
lFdsTEP.LT. 5) r,0 TO 74 
iFdSTtP.jrO.S) GO TO 152 
IF (MOD ( ISTE. p dO) .NE.O) GO TO 123 

15? tper=tpert=tr=tri=o. 

DO So I =2 . NP 1 

trp=tr 

TRsEmuT ( I) *0UDY (J ) *R (I ) 

TRI=TRl>.5tt(TR + TRP)«(Y(I)-Y(I-l) ) 

TPERPsTPER 
TPFRsRT W ( T ) <>TR 

56 TF'ERi=TPERI^.5-» (TPER^TPERP)*(Y(I)-Y(I-1.) ) 

Pl)MP=POE 

PDE=»TPERI/TRI 

IF (IsTEP.NE.5) PDE=.5« (PDE + PUMP) 
l'=2.f! 

ALPHAS. 55 

AH=Sin(3. 14l59»(PnE-.5) )-l . 

IK (PPE.LF. 1 . ) Al.PHA=. 55».213«AH 
G=1 .- ( 1 .- ALPHA opoE) /w 
0=0/ (1 . * (PDE-1 . ) /W) **2 
0=7.4 n7*o* ( 1 .-ALPHA ) /** 
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123 CMU=.oS<»g-.05:h«ACC 
r,o To 12l 

73 CMu=.oy-.n4<*ACc 

IF ( IpO.Eq.3) Cmus.09 
go To i2i 
7~» cmj=rMuiN 
121 RTCHu=SO«T (CMI.J) 

AKCOO=AK/rMUo»o.75 

uip=r<Ju,i) 

C *«*aXIAl VELOCITY GRADIENTS 
Dl.iDY c 1 > =0. 

DUDY/sgRl) =0. 

DUDYsO(l) = 0. 

DUDYsG(NPl) = o. 

VOrFp = Y (31 - Y(2) 

DUDYP = ARS ( (F ( JU , 3 ) - F ( JU, 2)1 /YDIFP ♦ l.E-30) . 
GO To <1,?,2), KIN 

1 DUDY (?) a (.O.S*YDIFP ♦ YI)/YI « DUDYP 
GO To 3 

? ouo y (?) a yi/<o.5<»ydt fp ♦' yu*dudyp , 

3 DUDYSO (2) = 0I>DY(3) # *2 
DO 9 l s 3» NMl 

YDJFiji ~ YDIFP ' : 

YOIFP =* Y V I ♦ 1 > - Y(I) 

DI.JDYM = DoDYP 

DUDYP a ARS ( (F ( Jtj, I* 1 ) - F ( JU . I )■) /YDIFP ♦ l.E-30) 
9 000 Y { I > * ( y ( I ) - Y < 1-1 ) ) / ( YOIFP/OUDYM^YOIFM/DUDYR) 
GO To (4,5,5) , KF.X 

4 OODY(M) a <0,5»YOIFP : ♦ YE ) /YE»DUDYP 

GO To 17 , • 

5 DODY(M) = YE/(0.5 < »YOIFP ♦ Y£)4DU0YP 
17 OODYSO(M) = UUO Y ( N ) *>**2 

DO 1? I =1 , NP 1 
OOOY.co ( I ) = QUOY ( i ) tn >2 
13 HTW(n=F(JO,I)/(F(jK, Ij ♦l.E-30) 

C <h»*TURO(JLENT VISCOSITY FORMULATION 
DO 1 0 3 0 I = 2»N ’ 

FUUT (T. ) =CV‘J tt HHO ( I ) «F( JK , I ) /HT'Y( I > 

1030 F (JUVf 1) =-E'-U)T (l) *fOUOY (1) /RHO(I) 

DO 1 033 1=2, NP l 

1033 TF(F(JU,I) .LT.F(JU,I-1) ) F ( JUV » I ) =-F ( JUV 1 1) 

C tx»«AD0 laMINaR AND TURBULENT COMPONENTS 
TF(ColT) 1050,1050,1052 

1050 DO 1 051 1=2, N 

1051 GAM< I )=EM(iLl I) +EMUT (I) 

GO To 1065 

105? DO 1 06 I=?,N 

IF(EmI)T (I)-EMUL (1)»CRIT) 1060,1060,1061 
1040 GAM M ) =Emi'L ( I ) 

FMUT'(T)=0. 

GO To 106 

1()61 GAM ( T ) =EM|)L ( I ) +EMUT ( I ) 

10/s CONTINUE 
106? m<T(l)eO. 

EHUT(MPI) o 0. 

C *««AXIM, PRESSURE GRADIENT (SOURCE TERM) . 

OAM3cGAM(3) 

GAMNbGAM(NMI) 



Q O O 


00 1 ; 1 Ial.NPl 
sU(i)=-DPox 

ill snu'isu. 

pftupn v , 

C+MMUMfTUpHi.iLENT KINETIC ENERGY 
299 iFU.NE.jt) <30 TO 399 
00 201 I=?.N 

?0l OAM ( I ) =EM|'iL ( I > /PM < JK) + EMUT (I) /PRT(JK) 

DO 2]? I = l.i'Pl 

on ( I) = F.mUT(I ><M.)UDYsO(I) 

21? SO ( I ) =-RTw ( I ) #RHO ( I ) 

IF (KJM.NE. 1 ) 00 TO 215 

FJK2=-KJI< JU) / (R ( i) orHO ( 1 ) «RTCMU) 

SO (2) si .EaO^FJK? 

SO(2)=-1.f30 
215 JF(KFX.NE.I) RETURN 

FJKNrRJE ( JU) / (R (NP1) «RHO (NP1 ) #HTCMU ) 

SU(N ) =1 .E3i.H»FjKN ■ 

SD(N >»~1.£30 
. PFTUPN 

C*<HHM>r)lsSlPATlnN KATE 

399 IF(J.nE.jD) GO TO <*99 
00 flnl Is?«N 

001 OAM ( I ) =K.MmL < I ) /PK ( JO) ♦EMUT ( I ) /PRT < JD) 

oo (Jin i=i,npi 

su ( I > saCt:i«EMUT ( I ) »HT»v ( I ) ( I ) 

CF2=1 ,92-,o667»ACC 
J F ( IPD • EQ . 2 ) CF.2 = 1.92-.1336*ACC 
IF { IP0.EQ.3) CE2= 1 . 9? 

MlO SD ( I ) =-CE?*KHO ( I ) «RT>i ( I ) 

IF (KTN .NE. 1 ) GO TO 812 
SU (2) =FJK-><H*1 ,5/AKCOO/Y (?) «1 ,E + 30 
Si) (2) = - 1 * E 3 0 
PI? IF(KFX.NE.I) RETURN 

SM(N )=FJKN <h»1 .5/AKCOG/ (Y (NP1 ) -Y (N ))<>!. E30 

SD ( N ) a -1 ,E30 

RETURN 

C***«#sTa0NATI0N ENTHALPY 
•499 IF ( J.mE. JHS) 00 TO 599 

<>#*PRANnTL/SrHMIDT MJMBF.RS FOR NEAR-WALL REGION (ROTTA) 

AX I SYMMETRIC FLOWS ONLY’ 

N.R. SAME RELATIONSHIP FOR SPECIES DIFFUSION 

IF (KraO.Eo.D GO TO 331 
IF(Kfy.NE.I) 00 TO 331 
00 330 1=1 ,NP1 

YW = 1 ,-R ( i ) /K (MPl ) 

330 SCH ( J) =.95-.4S#YW«YW 
00 TO 332 

.331 OO 333 1 = 1 * NP i 
333 SCH ( I ) =PRT (OHS) 

.33? DO 3fil I = j> * N 

301 OAM ( I) bKmmLI I) /PH ( oHS) +EMUT ( I) /SCH< I ) 

DO 3ri2 1 = 1, NP l 
.30? s D(I)«*F(JU»I )**2 
DO 303 I =? * Nm i 



FACU(SCh(I>-» . ) "(GAM ( I> ♦GmM{ I *1 ) ) / ( VDI F ( I ) + 2 , E-3 0 ) 

303 S D ( T ) =F Ar-l* (S 0(I + 1)-S 0<I))"0.5 
DO 3 1)7 I=p,Nf>ll 

FAC2= <SCH< I > /PRT ( JK)-1. )"<6AM ( I ) ♦ 6Ato { I ♦ 1 > ) / < YDIF ( I ) ♦ 1 . E-30) 
307 S 0 <n=s 0{I)*FAC2«(F(jk,I*l)-F(jK,I)) 

sD(l)aO. 

SD<N)=0. 

no 30/* IbjiN 

304 SU { I ) = ( SO *( I ) -SD ( I — 1 ) ) "2 • /TDA { I ) 

00 305 I = 1, NP 1 

305 SU(I)=0, 

SU(l)=i>. 

SU (Npl ) aQ. 

RETURN 

C"" 

C"""""MASS FRACTION OF HYDROGEN ELEMENT 
509 IF(J.NE',j‘a) GO TO 699 
DO 901 I = p ,N 

• 901 OAM ( I ) afMi'iL ( I ) /PR ( JA > ♦ EMtJT < I ) /SCH ( I ) 

DO 90? 1 = 1, NP1 
902 SU ( I ) =SD ( J > =0 . 

RETURN 

C""""" 

c """MEAN SUUaRE CONCENTRATION FLUCTUATIONS (OF H ELEMENT) 

699 IF( J.NE, Jr) GO TO 799 
yD I F p = Y ( 3 ) -Y (?) 

DCOYr = ARs'( (F (JA,3)-F(JA»2) ) /YDIFp+1 .E-30) 

JF (KlM.Ef). 1 ) GO TO 711 
DCDY = Yl/( .5 u Y0lFP-fYl)"DC0YP 
ncOYso=OCOY w "P 
SU(P) =CGl"DCl)YSQ"EMUT (2) 

SO (>») =-CGp"RTW <2) "RHO (2) 

SU ( 1 ) ='). 

SO ( 1 ) =-CGp"H r.l ( 1 ) "RHO < 1 ) ' 

NO TO 712 

7J1 5?= (Pjl (J) /H (2) /RHO <2) ) ""2"CG1 /CG2/F JK2 
SU (2) =G2"J ,E*30 

SO {?.) a-1 .6 + 30 
71? DO 7 J. 0 Is3»NMl 
YDIFmsYDIFP 
yD1FF=Y(I*1)-Y(I) 

OCOYf* = OCDYP 

OCOYp = AUS ( <F (JA, Itl) -F ( JA«I) ) / YD I FP ♦ 1 . E -39 ) 

DCDYSO= ( (y ( Ul)-Y (1-1 ) ) / { YU IFP/DCDYfU YD IFM/DCOYP ) ) ""2 
SU( I) =CG1 #dCUYS0"EmUT (I ) 

710 SO (I) a-CG2*RTW ( I ) "RHO ( I) 

IF(KFX.EQ.I) GO TO 713 

DC0Y=YE/(.5 # YDIFP+YE)"0CDYP 

DCDYsQ=DCDY*»P 

Sll(M) =CG1"DC0YS0"EMUT (N) 

SD (N) =-CG2 <>M TW(N) "RHO (N) 

SU (NP1 ) =0 . 

Sf)(NPl ) =-rG^"RTW (NP1 ) "RHO (NP1 ) 

RFTUpM 

713 C,M= (PJL ( J > /H (N) /RHO (N) ) "»2<KCG 1 / CG2/F JKN 
SU (N) =GN"l ,E* 30 
SO (N) *-l , f*30 
799 CONTINUE 
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s i 1 R r? p 1 1 t i im f stride 

DIMENSION AHL (4.1 ) , HOMT3 (40 ) * QMS (40) 'PliOM (40) ,RGQM(40 ) , THL <40 ) 
DIMENSION A (40) ,M (40) »C (40) ,1) (40) 

C0MI'iO('i/GEN'HAL/ACC*CSALFA.I)P()XflDX,ENTH (40) » F (10*40) , F S ( 1 0 , 4 0 > « GAM3 , 
1 f ; A M N . R A M ( /, 0 ) . I tlFINtlNDE ( 1 0 ) , INO I (T 0 ) , I STEP , IUTRAP , I TEST , KEX , K IN , 
?K-RAD,M»NE0.NP1 ,NM1 ,OM(40) » OMR (4 0) , BOM (40) » K AST ( 1 0 ) , PE I , PS I E , PS 1 1 , 

3 P (40) ,R(4o) » HHO (40) » R J E ( 1 0 ) ,PJI (10) , KME , KM I , SO (40 ) ,SU<40) * WM (10) . 
4XD,X|I, Y (40 ) » Yi'IF < 4 0 ) * YE»YI ,RJTE (10) ,RJTI (10) 

COMMON/Cjc/JU, JK, Jl)» JHS, JA, JG, JTE. JLE . JUV* JB 
POMMoN/CFs/JHP. JOP, JOH. JH?0, JH» JO» JN2 
COMMON/AUXST/ TOa ( 40) ,U1P 
POMMoM/WFcT/FDlFS (20,2) ,TS(20*2),BP(2) 

COMMoN/CPR0P/lO2, IN2, I H20 » OFAC » NR » NS , WQO , WRO, WSO» WTO , GASCON, GAMMA 
STRIPE (0) 

C „ <H»«pONTRnL INDICES 
ENTRY BTRtDO 
NPl=M+i 
NM1 = N-l 
P M ( 1 ) n 0 , 

OM(NRl) =1. 

G A M 3 .. = 0. 
gamn = o • 

ISTEf' = 0 
TFIN-o 
T U T R A P = 2 

c <h»»zfroimg of Important arrays 

DO 3 f>3 Jsl»N£0 
kast < j) = n 

353 RJE ( J) *HJT ( J) =0. 

00 is 4 1=1, MR! 

DO 3^4 J = 1 , N S 

354 FS(J,I)=0. 

DO 355 1=1 » NR 1 
DO. 355 J=1 , 10 

355 F(J,T)=0. 

DO 350 1 = 1, NR 1 
354 F N T H < T ) = 0 . 

.RETURN 

C«hhh> sTRInF (1) 
e owttOMEGA RELATIONSHIPS 
ENTRY S T R y D 1 
DO lOO I = 1. N 
Ol'D ( T ) «0M( 1*1 ) -OR ( I ) 

100 OHS ( T ) =0M > 1 + I ) ♦OR ( I ) 

[)0 10) I = 3. NMl 
BOM ( T ) =OM ( I+I ) -OP ( 1-1 ) 

ROMT3(I>=BOM(Im>3. 

)01 PONT I MU£ 

HOM(p) = pM ( 3 ) + OM ( 2 ) 

HOM(r-i) = p. - OM(N) - OM(NMl) 
nop = OM ( 3 ) /Of'O ( 2 ) 
n r > 3 = 1 . - omr 

01 'S 2 a OMj3)<m»p 

OMS3 = Om (?)**? 

OMS2 c OMSH/ (OMS2-OMS3) 

0MS3 = 1, - OMS2 
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OMN = <OMf)(N) ♦ OMD (NM1 ) ) **2 

OKNMi = Of/0(N)«<»2 

OMN = OMN/ (OMN - OMNMl) ‘ 

OMNMl = 1. - OMN 
Y ( 1 ) = 0 . 

BP(l')el. 

RP ( r! j =1 •' 

JNDE { JU) =1 
I HD I ( JO > = 1 
RETURN 

STRIPE (?) 

PNTMY STR j 02 

C *»«TEST for negative velocities 
101 o IF (IiiTRAP.EQ.O) GO TP 1011 
00 1012 I a 2, N 
TF(F(JU,d .GT.O.) GO To 1012 
WRITF (6, 1?00) F(JU,I) ,1,1 STEP 

1?00 FORMAT (/1H0** NEGATIVE VELOCITY 0F**lPeil.3»» CALCULATED AT NOD 

1E#»I3.* AT 1STEP=«,I5) 

F ( JU, T ) =1 .E-30 
IF(IUTRAP.GT.I) IFIN=1 
101? IF (IllTRAp.uT.2) ITESTsl 
C »<H> c Ro S S _ 3 T R p A M DISTANCES (Y*S ANO RiiS) 

C *##P0GE REGIONS 
1011 RHP s RH0’(2) * F(JU,2> 

RUlaRHO ( 1 ) <U r ( JlJ. 1 ) 

RURAJ a Rl'll/WUP 
00 TO (1013*101/+, lolA) , KIN 
10 14 GO Tn (101fl*l019)» KM AO 
101 a B P ( l ) = 0.333333 ♦ 0 . 666667 »RURAt 
00 TO 1013 

K-iQ RP(l) a (R ( 1 ) # (5.<»HURAT+1 . ) ♦3. tt R (2) * (RURAT + 1 . ) ) / (R < 1 > *R (2) ) /6. 

1013 Y I = pEI a OMP ( 1 ) / (UP ( 1 ) «*RUP) 

RUNP] = RHO(NPl) « F(JU*NP1> 

RUN c RHOfN) * F ( JU * N ) 

RURAT c RilNPl/RlJN 
GOTO (1'0?0*1 021*1021), KtX 
1021 00 To (1024*1025), KRAD 

1024 HP { 2 ) = 0.333333 + 0 . 6666b7«RURAT 
GO To 1 02n 

1025 P P ( 2 ) a (R (NP1)*» (5.*RURAT + 1, ) +3.«R (N) » (RURAT+1 . ) ) / (R(NP1) *R (M) > /6. 
1020 YF = PFl*oMP(N)/(RP(2).#RUN) 

<MH>y /s , r/S, TUAZS ANO YDIF*S 

Y0IF(I)=2. W (Y(I*1)-Y(I))/R (J+1.5) 

Y (2) = Yl 

YDIF(l) = 3» * YI 
PO 1 0 1 T I = 2, NM1 
TDA(i)=PET*60H(I)/RUP 
WUMaptlP 

RUP = RH() ( 1 + 1 ) *F ( JU, 1 + 1 ) 

Y I > T F ( I ) ap r I*OMn <!)<>( l./RUM + l./RUP) 

1017 Y(1.+ 1)=Y(t)+YUIF(I)». 5 
TP A (r)=PET«UOM ( N ) /HUP 

Y ( N P 1 ) a Y (N) ♦ YE 
YDIF(N) = 2. « YE 
IF (KPAD.Fc 1 ) RETURN 

r. «#*Mon if I gaj t ons for axial symmetry 
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11?3 tF{CsaLFa.FU.O.) GO TO 1110 
C «**CSALFA .NF. ZERO 
COSO?* • 5#rS^LF’A 
TF(W(1 >.Nf»°*> GO TO 1105 
c <n»#R(l) .EQ. ZERO 

DO UOb I = ?-■> NP1 . .. 

Y<I>sS'WT<AHS(Y(I)/GOSD2)> 

110a R ( I ) cY U ) «CSALFA . 

00 TO 1107 

C <h>»R ( 1 ) .RE. ZERO 
1105 P1D2=.5<>R', 1) 

RlD2sO=WlD2 1> Rin2 

DO 1104 I = ?-* NP1 . . . ‘ ' . ' . 

Y ( I ) =Y (I ) / (R 10?+ SORT (A8S(R102SQ+C0SD2*Y (I) ) ) ) : , . 

1104 R(I)aP(l)*Y<I)*CSALFA 
1107 DO HOB I = 1. N 

llOn YDlF(I>«Y0lF(I)*4./(«(i)+R(I*l) >##2 ■ . ■ . • . 

GO TO Hi? 

C. •ottttf.SALFA • Ep • ZERO . . ; 

1110 DO 1)11 I = ?.i NP1 
Y ( I ) =Y < I ) /R ( 1 ) 

1111 RllbPlll : . ' 

R1S«=P<1>«* 2 • ' . , . " 

DO 1 1 0‘> I=2iN 
nog YDIF(I)=YDIF(I)/R1SQ 
11 1? YI = Y (2) 

vf: = Y (NPl ) - Y(M 

RF'TUPM ... 

C **#« STRIDE (3) 
c <hmvmaiN NUIIeRiCAL method 

FNTRY STrtD'3 , - ■ 

IF (KIN . r(3 • 1) CALL WF<1) 

IF (KFX ,FQ. 1) CALL Vi F (2) 

C ■»<K»PRELT , -ll N Af ? I £ S 
3000 GsRMj-UHE 
P X = P E I / 1.) X 
PD4 = , ?5 h»PX 
r,t)Aa,PP*G 
PGDAePDA^oOA 
RMlD?=«b<»RMl 
PG0M(1)=0. 

c AHL • PGOM AND RBOM 

THL<1) =0. 

DO 3010 I a 2. N 
HL=RmTD2-gD4»0ms (I) 

THLU) = 2.”HL 
AML (I) = A05(HL) 

PGOM'(I)=PGD4-#Onnm 
3nlO.PR0*(I>=PX»«Ot4m 
C •»#«STaRT OF j loop 
DO 33?0 J=1»'NE0 

I F ( J , ME . J 1 1 ) CALL AUX(J) 

TTP=0. 

PGOHP a 0. 

OAM(NPl) s 0. . 

c <»««settjno up coefficients 

AKAST s FLOAT (1 -KaST(J)) 
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KASTl = KiSS'f (I) ♦ 1 
INDEX I a lNDI ( J) 

indexf=Indf < J) 

DO 3004 I = 2 » N 

su 1 1 ) = sm(I) <* TDA(I) 

SD ( I ) a SD(I) * TDa<1) 

PGOMM = PgOMP 

PGOMP = PoOMU) * AKAST 

TTM=TTP 

TP = ((5AM (I) ♦ QAM (1*1) )/YOIF(I) 

TTP = TP * AHL(I) ♦ ABS (TP-AHU I ) > 
a ( I ) =TTP-THL 1 1 ) -PGOMP 
B (I ) sTTM+THL (1-1 ) -PGOMM 
T ND 1=?/I 
IMD2=1/(NP1-I) 
lNDtX=l+INDl*IND2 
GO TO (303*30?) , INDEX 
303 GO To (304*305), KASTl 

30a c ( I > =PD4«'(B0MT3 ( I ) «F ( J« I ) ♦OMQ < I ) «F ( J* I* 1 > *OMD <I-1)*F(J»I-1)-)*SU(I) 
GO To 306 

305 C(I')cPtfOMf I>*F(J*I)*SU<I> 

306 D ( I > sPOOM f I ) “SD ( I ) 

C #*4moDIFICATTON 3 FOP BOUNDARIES 
jMDEX=l+INDl*2»lNn2 
GO To (3004»300B»3002) ♦ INDEX 
3 0 Ofl A (I)cA(I) +PGOM(2)*AKAST 
B ( I ) =?.«RmI 

GO TO (30n9»3"04,3004) , KIN 
3009 GO To (3011*3003), INDEXI 
3011 TT=2.»TS(J,1) 

B ( DsAMAXl (TT*B (I) *0.) 

C(I)eC(I)-TT*FDlFS(J,l) 

GO To 3004 

3003 B ( I) =0. 

D ( I > sf) ( I > *2 • I 
C(I)nC(I)*2»*RJTI(J) 

GO To 300 a 

300? B ( I ) sH ( I ) +PG0M{NP1)<* AKAST 
A ( I ) s-2.oHME 
GO To (3012*3004)* «EX 
301 ? GO To (3013*3005)* IND£XE 
3013 TT a 2.»rS(J,2) 

A (I)aAMAXl (TT + A (I) ,0.) 
rm=c<i)-TT#FniFS(j,2) 

GO TO 3004 
3005 A (I)so. 

D(I)sD(I)-2.*PM£ 

C(I)»C(I)-2.*MJTF (J) 

3004 D ( I) aO( I) + A ( I ) +B ( I) 

C 4«*AOJUc;T FrfE-BOUNDAPY VALUES 
Rill a RHO ( 1 ) » F(JU*1) 

IF (KJM. Me. P. OR. RU1 .E/J.o. ) GO TO 3006 
F ( J * 1 ) a (F(J*l)*SUU)*nX/RUl)/(l.-SD(l)*OX/RUl) 

30 06 RUNP]sKH0'<NPl) «F ( JlJ»NPl) 

IF (KFX.NF..2.OR.RUNP1.E0.0.) GO TO 3007 

F ( J»MPl > a (F ( J,NP1 ) *SU (NP1 ) <*()X/RUNP1 ) / ( I . -SD ( NP 1 > "0X/KUNP1 > 

C ■»<M»90LVp FOR DOWNSTREAM F*S 
30 07 B(2)rs{B(2)*F ( J,l) *C (?) )/D (2) 


66 


A (2) =A (2) /L>(2) 

DO 3o2l I a 3. N 
T=0( p-B(T)*A(I-l) 

A ( I ) sA ( I ) /T 

•3021 0II)s(B(I)«H(I-1)*c(I)I/T 
DO 3o?2 i o ASHa 1 , N'M 1 
T=NP1 -IOa^H 

302? F(J»l)=A(n«F(J,I*l) >B(I) 

c **#adjijst boundary values 

GO TO (3210*3220,3230), KIN 

3210 GO TO (300,3211 ) r iNDEXI 

3211 RMT=PMl*To ( J * 1 ) 

RJI(J) = RJTI(J) - RMl#F(Jtl) 

IF (PUT .pT. 0.) GO TO 307 
F(J»1) = F (J,2) 

GO To 3226 

307 F ( J, 1 ) = (RJTI(J) ♦ TS(J,l)*(F(J,2) + FDIFS ( J* 1 ) ) ) /RMT 
GO To 3226 

300 RJI(J) => TS ( J » 1 ) # <F (J,1)-F(J,2)-FDIFS(J»1) ) 

RUT I ( J) a RJI(J) ♦ RMI«F(J,1) 

GO To 322n 

3230 IF(R(1> .Fb.O.) GO TO 309 

F { J «"1 )=F (J,2)<»0MS2+E(J*3)#0MS3 
GO To 322n 

309 F(J»i)aF(J,2)«>OM2 + F(J*3)«0M3 
3220 GO TO (3310,3320*3330) ♦ KEX 

3310 GO To (3ln,3311) . iNOExE 

3311 RMT=-RMF+TS <J*2) 

RJE(J) = RJTE(J) - RME«F(J,NP1) 

TF (RMT ,pT. 0.) GO TO 317 
F|J,NPl) r F(J,N) 

GO To 332n 

317 F(J,NP1> a (-RJTE(J) ♦ TS (J,2) <» <F < J,N> ♦FDIFS < J*2) > ) /RMT 
GO To 332n 

310 RJE(J) ■ -TS(J,2)tt(F(J,NPl)-F(J,N)-FDIFS(J,2) ) 

RJTF ( J) = RJE(J) + RME*»F ( J , NP 1 ) 

GO To 33 2o 

3330 F ( J.NP1 ) sF ( J»M)#0MN*F ( J,NM1)*0MNM1 
3320 CONTINUE 
XOaXO 

ACC»'(F(Ju’.1)-U1P)/DX 

PSII=PSII-RMI»OX 

PSIFaPSIE-RME^DX 

PEI=PSIE-PSII 

ISTEPaISTFP+1 

return 

END 


COMPTLirR SPACE 


SU8H0i.lT INF '"'F(k) i 

C0MM 0 N/GFNRAL/ACC,CSaLFA, 0PDX,DX,ENTH<40) tF(10.4fl)*FS<10.40) *GAM3, 
1GAMN,6AM( 4 o) »I,IFIN» lN0Etl0> ,INDI(10) . I STEP , IUTRAP , 1 TEST , KEX , K IN, 
?KRA0.N*NE0»W1 »NUl.OM<40> • OMD (40) tBOM(40) »KaST(10) »P£I *PSIE»PSI 1 * 
3P (40 ) »R (4n) » H H 0 (40) .HJf(lO) * R J I (10) ♦ RME ♦ Rrll , SO ( 40 ) ,SU (40) ,WM ( 10) » 
4XD,Xu. Y (4n) ♦ YU IF (40) . YE«YI,RJTE < 10 ) ,RJTI (10) 
rOMMOM/CJc/JU. JK, JO* JHS* JA« JG» JT£* JLE» JUV» JB 
COMMOM/CF c/OHP , J02 , JOH , JK20 , JH , JO » J NP. 

r0MM0N/CTr-H8/AK,ALMG.CMtl,qMUIiM»CEl *CE2tCGl » CG2 , CR I T * DUDY < 40 > » 
lDUOYsQUo) ,EMUL (40) »EMUT(40) » IPO, PR (10) , PRT (1 0 ) * RTW ( 40 ) 
rOMMOM/CWwF/YH (?) , IPUF (2) *EWALL*H 
COMMoM/WFsT/FOJFS(20,2) »TS(20,2) ,8P(2) 

OOMMON/CI^MA/SCH (40 ) 

C ###########»»###« »«»«# 4*# »»!>'»»»»« » 

DATA SUALF/-04/ 

00 To UOi U ) »K 

1 0 I'-'=l 

1 Ns2 

GO TO 12 

11 J U'aNP 1 
I N=N 

J? CONTINUE 

C 4#*RFFEPFNCE QUANT I T ies 

|]HFF = A8S (F (UU, IW) -F (JU»IN) ) 

RhORFF=O.E»(RHO(IW)+RHO(lN) ) 

RUREFsRHOREF<HJRFF 

RPEFaR(lN') 
krurff=rrfF # puref 
viTeFscmul '( f i w ) 

YREFaAHS ( v ( lW) -y ( IN) ) 

RF = HuREF#yREF/VRF.F 
AKa (PM I - (RME -RHI ) <*OM ( IW) ) /RRUREF 
FFs-SU ( I W ) «yrff/rijref/uref 
12? TF (RF.LT. 132. 25) 60 TO 110 
c o«<m_og law assumption 
121 CONTINUE 
LAM = 0 
NIT = n 

101 SHAI_FI=SHaLF 
yrp=rvk<>shalf 

c «*4CALCul.ATloN OF e FOP ROUGHNESS 
F=EWALL 

I F ( I pnF ( K ) . EQ. 0 ) Go TO 16 
c 4 **sand-grain Roughness 

IF (YRP.GT.3.3333) E=30./YRP 
If, FRaRpoE 

S = SH/\LF*4? 

SLOCaS ♦ AM+EF 

IF (S|. OC.GT.O.) GO TO 104 
SLOCal .F.-30 
SHALFaSORT (ABS ( AM*EF) ) 

104 BEE=SORT (sLOC/AK) 

AR(J=FW* (SHALF* (AM/ (1 .♦BEE) ♦,5«EF) /SHALF) 

JF ( APG.GT. 1 1 .<>£) GO TO 106 
r,0 TO 110 

104 SHALFbAK/aLOG ( AHG) 

JF(ApS<SHaLF-SHALF1) .LT..0001.0K.NIT.GT.10) GO TO 102 
NITbMIT*! 
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GO TO 101 
lo? s=shale**5 

BP(K)=l./' f l. ♦)}££) 

GO TO 103 

C *#*LAMIkaH FLOW 
lln LAM = 1 

AMREc AM*Rp 

FRE=EF«RE 

I F ( A G S ( A MR E ) * L T . , 0 1 ) GO TO 111 
FXPMrF = FXP ( AMRfT ) 

STORFqEXPmRE-1 ."-amre 
AMRESO=AMRE*AMRf 

SRE = AMRE»'c 1 ••?STOWE*F'ME/AMHES(J) / (E.XPMRE-i , ) 
0UTlrS«F>sT0HE/AMRESQ*FHE« (STORE-. 5*Ai*.RESQ)/(AMRES0»AMRE) 
GO TO 112 

111 SR6=(?.-FRE <> (l.+AMRE/3.))/(2.*AMRE> . 

OUTlrS«E«(f.5 + AMRE/6.) +FRE*(.16667*AMRE/24.' 

11? TF (SPE.6T1.E~30) GO To 113 
SPE=].E-3n 
0UT1=. 33333 
1.13 .<? = S«F/»'<E 
_ ‘ RP ( K ) aOUT 1 

103 00 SOOO j= 1 i NEQ - 

EOIFs(J,K)=0. 

IE(J.FU.JHS) F0IFS(J.K>=(H-1,)«.5«UREF*»2 
TFcJ.mE.Jh) GO TC 200 
TS ( J.K) =SwRHUREF 
GO To 500n 

c «m»#stagmaTion Enthalpy and mass fraction 

POO CONTINUE 

IE (Rr.LT, 132. ?5) GO TO 210 
IF (L/.M.EO.l) GO TO 210 ■ 

POl CONTINUE 

PRRAT = PW < J) /PUT (j) 

IF (KpaO.Ep.2) PRRATsPR (J) /SCH (N) 

PJAY=:9.^ (PRRAT-1 . ) /PRRAT»<*.25 
SF = S/(PRT'(J)« (l.*PJAY«SHALF) ) 

JF (KraU.Eo.2) SF=S/(SCH(N)»(1.*PJAY#SHALF) ) 

GO To ^13 

C «**LAMlMAR flow 

Pin TF(Aps(AHRE> .LT..01) GO TO 211 
SF=AM/(EXP<PR(J)»AMRE>-1.) 

GO To 212 

211 SF = 1 ,/PW (J)/RE/(l .*.5»PR (J)*AMRE) 

PI? CONTINUE 
? 1 3 TS (J.K) sSF^RRUREF 
SOOfl CONTINUE 
RETURN 
END 

comptlfr space 
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SUBHni'T INF LAMGLl : V 

DIMENSION CHMW(7,60) ,PTfc'MP(40j ,RCON(10) » STORE (10)* HEM (10) 
(-OMMqM/GfNRAL/AcC.CSALFA ,DPDX.OX,E^ rH(40) ,F (10,40) ,FS( 10.40) ,GArt3, 
1 rAMN , GAM ( 4 O) ♦ I. IFIN, INDE (10) ,INUI (10) . I STEP , I UTRAP , I TEST , KEX ♦ KI h , 
Fl<RAO,N. Meo.NPI »NM1 , Of' (4 0) »()MO (40 ) .ROM (40 ) * KAST ( 1 0 ) . PEI » PS If . PS I I . 
3P(40) ,rt(4n> * RH() (40) .PJEO'M ,RJI (10) . RME ♦ RrtI . SO ( 40 ) ,SU(40) ,WM(10) , 
4XD,X(i, Y (4n) . VO IF (40) , YE » Y I « RJTE (10) ,RJTI (10) 

C.OMMON/Cjq/JlJ, JK. JD, JHS, JA, JG, JTK.JLE, JUV, JB 
rOMMnw/CFc/JHP, jn?, joh, jh20, jh, jo, JN2 

COMMOM/CST AH/MO (7) , CPBAR (7,60) »HT (7,60) .RC (5,60) ♦ ISTAR (40) , F A CM 
COMM0N/CTmHB/AK,/)LMG,CPU,CMUIN,CE 1*CE2.CGI , C<32 , CR I T » DUDY (40) . 
lDHDYjO (40) ,tMUL (40) .EMUT (40) . IPO, PR ( 10) ,PRT (10) ,RTW (40) 
rOMMnN/CPR0P/I()2, I NR, I H20 , OFAC . NR , NS , WQO.WRO.WSO, WTO , GASCON , GAMMA 
C iHKt # # <MM» tt # 

CHAPTER 1 

C«<MHHtLOA(3iwu OF EQUILIBRIUM CONSTANTS, ENTHALPIES AND SPECIFIC HEATS 
p N T R Y LANrLI 

C *<M>RfcFEPFNCE ENTHAI.PIES 

DATA (HO(J) »J=1 ,7) /-2023. 8.-2074.7,72 05.6, -60164.7,5061 6.5, 57949. 1 
1,-207?. 3/ 

C «•»«(- H £ M j C 4 L pOOILIRRIuH CONSTANTS 
DO 2(1 J = r,NR 

READ ( 5 , lOnO) (RC(J,I> .1*1.60) 

DO 20 I = r,60 

RC(J.T)=APIN1 (RC(J,I> .100.) 

?o RC(J,I)=lo. tt<> RC(J,I) 

C ###C0NSTANTS FOH GLOBAL REACTION 
DO 15 I = l'. 6 U 

15 RC(5,I)aRr (4,I)*PC(3,I)/C«C(1,I)*S0RT(RC(2.I) > > 
c <MH» C 0NV F RSI0N INTO CONSTANTS FOR CONCENTRATION RATIOS 
A10 =VM( JHp)/WM(JH)#»? 

A?0*WM( JO?) /^M(JO)**P 
A30=W^ (JOH) /WM ( JO) /WM ( JH) 

A4 0 = (-.'M ( JHfO) /PM ( JOH) /WM < Jn) 

AbO=HM < JH?0) /WM ( JH?) /SORT (WM ( J02) ) 

DO 571 Is'l.bO 
RC ('l.r)eWo(l* U«A 10 
RC (2, T>=R 0 (2,I)»A20 
RC (3. I ) =Rp(3, I ) «430 . 

PC (4, T ) =Rp ('♦* I )«A40 
571 RC(5,I)=Rb(^>*I) tt A50 
C *«*FNTHALPIF.c AND MEAN SPECIFIC HEATS 
DO 3o J= 1 '.NS 
CFAs4)U7./WM(J) 

HO ( J) =H() ( j) #CF A 
READ(5.20n0) (HT(J, I) .1=1.60) 

DO 30 I = l',60 
HT ( J* I ) =HT (0,1 ) »CFA 

30 CP8AR(J,I)=.01#(HT(J,I)-HO(J) )/FLOAT(I) 

«<n » iteration parameters 

IPEG1N=3 
TTMAX=12 

cc=.oi 
RPC=,5 
RPD=l.-RPr 
C 
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FST0iCH=2.»wm< JH?) /'wm < J02) 

NRPsf.iR* 1 

HFTURN 

CHAPTFR ? 

C**#** INITIAL FNTHALPIES 
Fl'JTWY LAN«L2 
FNT H ( I ) =0 . 

FACX=.01#F(w)TE,I) 

IFaCX=Facx 

FACS=FACX-FLOAT (IFACX) 

DO 1 15 Jal ,NS 

FNTLcHT (J - , IF ACX) ♦FACS^ <HT (J, IFACX* 1 ) -HT (J» IFACX) ) 

115 FNTH(I)sENTH(I) *FS(JiI)«ENTL 
RETURN. 

C********** ********************************************. ***************** 
CHAPTER 3 

C*****LOCAL mass FRACTIONS 
ENTRY LANpl.3 

c ***locatf position where hydrogen concentration is stoichiometric 
X=1.-F(Ja',1> 

X=X*oFAC 

IF(F(JA»1)/<X*1.E-10) .LT.FSTOICH) GO TO 625 

DO 6p] 1=1, NP1 

X=I.-F(Ja'.I) 

X=X*pFAC 

PARA=RAT-1 ,E-10 

RAT = F(UA,if)/(X*l.E-10) 

IF (Rat. LT.FSTOICH) GO TO 1060 
API CONTINUE 
1060 TLOCp=I-l 
ILOCr=I 

FACMsY <I-1) + (Y(I)-Y(I-1) )/ (RAT-PARA).* (FSTO I CH-PAR A) 
go To 626 
625 Tl.oCpxO 
I L0CP= 1 

6?6 DO 2f>fl JOm=1»NP1 

t=npi*i-iom 

TF(I.LT.ILOCP) GO to 740 
GO To 74 2 
740 I=IL0CP-I 

74 ? TSM = AMAXl',F <JTF,I) ,250.) 

C 4**CONCFMTRaTIONS OF OXYGEN ELEMENT AND NITROGEN 
XN2Pn=l.-F(JA,I) 

x=oFac*xn?po 

FS ( JM?» I ) =XN2P0-X 

RAT=F (UA,t)/(X+1,E-10) 

DO 304 L = 1 , N S 
304 5TORF(L)=o. 

C ***FOUlL TBR I mM CONSTANTS AND MEaN SPECIFIC HEATS FOR UPSTREAM STATE 
FACX=.0l4TSM 

ifacx=facx 

FACS=FACX-FLOAT( IFaCX) 
no loss l=i*nrp 

1055 RCON (L) *Rr (L, IFACX) +FACS* ( RC (L » I F ACX* 1 ) -RC (L , IFACX ) ) 

CF AC = oHo ( T ) *GASCON*F < JTt » I ) *1 .E-05 

DO 1056 l=i»nr 

1056 RC0N(L)=RC0 N <E)*CFAC 
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WTON (MHP) sRCON (NRP) <>sQRT fCFAC) 

00 3) 7 Lai, NS 

317 cPmN(L» D.aC^BAR (L» iFaCX) *FACS# (ChBAR (L, IFACX* 1 ) -CPBAR (L. IFaCX) ) 
TTFH=rt 

IFf.c'OT. »T02.t;o.»).ANO.IH20.EO.O) ) GO TO 175 
C ^SOLUTION oF QUADRATIC equation For HYDROGEN-NITROGEN mixtures 

FS{JH»D=SQWT<1.44.*r<JA.I)*RC0N<l) j 

FS{JM,I)=.B tt (FS(JH t I)-l.)/RCON(l) 

FS<Jh?»I)=F.<Ja, I)-FS(JH, I ) 

00 Tn 208 . 

175 CONTINUE 

C *<*«STaRT up iteration cycle 

C ***DF.TEPMINE LARGEST TERM IN EACH GROUP 

c CQUArrUN FOR OXYGEN ELEMENT 

J 1 = Jo? . 

Ff-'1=F5U0?,I) 

lF(Fs(UO,r) ,LT.FS(J02,I) ) GO TO 210 
J1 = JO 

FM1=F5( J0'. I) , 

•?lo IF(FS(JH;» n ,I)«wsO.LT.FMl) GO TO 211 
J1=JH?0 

FH1=FS < JH?0» I H>NSO 

211 IF <FS (JOH‘, I)*w*ro.LE.FMI) GO TO 350 
J1=J0H 

FMlsFS < JOH, I ) «WR0 
35n CONTINUE 

C —EQUATION for HYDROGEN elfment 

J?=JH? 

FHPaFS UM?»I ) 

jF(FS(JH,f)«LT.FS(JH2,D) GO TO 212 
J? = JH 

FU2=FS < JH'. I ) 

?12 IF(Fs(JH2c,I)«W00.LT.FM2) GO TO 253 

J2=JH?0 

Fii2 = FS ( JH?0* I ) *WQO 

253 IF {FS(U0H'.I)*wT0.LE.FM2) GO TO 353 
J?=J0M 

FH2»FS < JOH, I ) #WTO 
353 CONTINUE 

PRECa 11 r Ion IF H?c IS Largest OF BOTH GROUPS 
TF { .NOT. (J1 .EQ. JH20.AND.J2.EQ.JH20) ) GO TO 279 
IF (HAT*G£.FST0ICH) Go TO 356 ! 

J 1 = Jo? : j- 

FM 1 =F5 ( JO? » I > 

IF<F$(J0,t).LT.F5(JO2,I)) GO TO 622 
J1=J0 

FM1=FS<U0'. I) 

H22 JF (Fs (UOH - . I ) # R«0,LT.FM1 ) GO TO 379 
JlsJOH 

FM1 = FS(J0H,I)<»WRC 
GO To 3 79 
35A J?= JH? 

F H2=FS ( JH? , I ) 

JF (F? ( JH, f) .LT.FS (JH2* I ) ) GO TO 361 
J 2=JH 

FM2-FS ( JH, I ) 

361 IF<FS(U0H,I)«WT0,LT,FM2) 'GO TO 379- 
J2aJ0H 
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FM2-FS < JOH , 1 > *WTO 
GO To 3Tg 
2"?9 CONTINUE 

C«»* PRECAUTION' IF OH IS LARGEST OF ROTH GROUPS 

IF (.COT. (Jl.EO. JOH.ANO.JP. EG.JiiiJo) ) GO TO 379 
IF (RAT.GE.FSTOICH) GO TO 367 •- 

Jl =Jr? 

FM1=FS< JO?,D 

IF(Fs<J0,t).LT.FS(J 02.I»> GO TO 358 
Jl=Jo 

FMIbFS < JO. I ) 

35* IF (Fs(JH2n»I)*WS0.LT.FMn GO TO 379 

J1=JH?0 

FNlaFS (JH20*I ) *tvSO 
GO To 379 
367 J?=JH? 

FwpaFS IJH9, I) 

IF(Fs<JH,t).LT.FS(JH2.I)) 60 TO 381 
J? = JH 

FU?=FS ( JH'. I ) 

381 IF (Fs ( I ) «WQ 0 .LT.Fm 2 ) GO TO- 37.9 
J?=J H ?0 

FM?=FS ( JHpO * I ) *W00 
379 rONTjMUE 

C «hh»oVFR-RIDING PRECAUTION' 

IF (F ( JA, I ) .GT. .9899) J1=JH20 
IF(F(JA,I) .LT.. 0001) J2 = JH2(V . , 

C — - — - HP OR : H AS J2 

IF ( J?.Nc. JH2.AN0. J?.NE. JH) GO TO 240 
IF (JP.El.JH2) F $ ( JH* I ) =SORT (FS ( JH2 » I ) /RCON (1) ) 

IF (Up.tO. JH) F S ( JH2 , I ) =FS ( JH* I ) ##2#RC0N ( 1 .) ■ 

IF ( Jl .to. JH20) GO TO 233 
IF (J) .EG). JOH) GO TO 427 

IF ( Jl . r-G , jo ) FS ( J02* 1 ) =FS ( JO, I ) »*2<>KC0M2) 

IF ( J] .LQ.J02) FS < JO* 1 > =S«RT (FS ( J02, 1 ) /RCON (2) ) 

FS ( Jh?0» I ) =FS ( JH2, I )*SGRT (FS ( J02, I ) ) «RC0N(5) • 

FS ( JOH* I ) =FS ( JO, I ) *FS ( JH. I ) <>RCON (3) 

GO To 890 

233 FS(J0?*I)=(FS(JH20,I)/FS(JH2,I)/HC0N(5) )«o2 
FS(JO, 1 > = C Q H T (FS ( J02.I) /RCON (2) ) 

FS (JOH* I ) =FS ( JO, I ) *FS ( JH, I ) *RCON (3) 

GO To 890 

427 FS ( JO. I) =FS (JOH, I )/FS(JH. I) /RCON (3) 
FS(J0?*I)=FS(J0,I)**2 # RCON(2) 
r S ( JH2 ( 'J, I ) =FS( JH, I ) *FS (JOH, I ) <>RC0N (4) 

00 TO 290 

?4n IF ( J?.N£, JH20) GO TO 250 

C — — h20 AS J2 

IF (Jl. tO. JO) FS ( J02 * I ) = FS ( JO , I ) *#2#RC0N (2) 

IF ( Jl .EQ.J08) FS ( JO. I) =SORT (FS( J02, 1 ) /RCON (2) ) 

IF ( J] .NE. JOH) GO TO 430 

FS ( JO, I ) =FS ( JOH, T ) * * 2 / F S (JH20* I ) «RCOIM (4)/ RCON (3) 
FS ( JO?* I ) =FS ( JO, I ) »<*2«RC0N (2) 

430 CONTINUE 

FS ( JU? ♦ I ) =FS (JH20, I ) /SORT (FS ( J02, I ) ) /RCON (5) 

FS ( JH. 1 ) =cQKT (FS ( JH2, I) /RCON (1 ) ) 

IF ( Jl .to. JOH) GO TO 290 

FS (JOH* I ) =FS (JO. I )*FS(JH. I)#RC0N(3) 



c 


r ,0 TO <*RO 


OH AS J2 . 


?50 cONT j NU£ 

TF (Jl.tQ. JO) FS (JO? 1 1 ) =FS ( JO, I ) «»2#RCON (2) 

IF(J] ,Eu.jo2) FS ( JO* I ) sSORT (FS (J02t I ) /RQON 12) ) 

IF < J\ ,KO. JH20) GO TO 501 

FS(JH. I ) =FS < JOH, t )/FS ( JO. I ) /HCON (3) 

F S ( JhP * I ) =FS ( JH, I )«*2<*RCON ( l ) 

FS < JH?0* I ) =FS ( JH. I) #FS ( JOH, I ) *RC0N(4) 

00 TO 2 l *0 

SOI FS ( JH • I ) sFS ( JHPO • I ) /FS { JOH « I ) /RCON (4) 

F? ( JH2* I ) =FS < JH, I )*h»?«RCON ( 1 ) 

504 FS(JO*I)=FS(JOH,I)/FS(JH,n/RC0N(3> 

FS { Jo?* I ) =FS ( JO, I )»«2»rtC0N (2) 

?9 0 CONTINUE 

C «hm»RE-EVALUaTE LARGEST TERMS FROM ADDITIVE EQUATIONS 
IF (WaT.LE.FSTOICH) GO TO H10 
030 CONTINUE 
Fm1 = x 

TF ( J1-.N2. JOH) Fm1=FM1-Fs ( JOH, I ) »WRO 
IF(Jl .HE, J0«) FM.lsFMl-FS(J02, 1) 

IF(JI.^E.jo) FM1=Fm1-FS(J0,I) 

IF { J ] .NF..JH2U) FH1=Fm 1-FS(JH20»I)«»WS0 
FS(Jj.I)sFMl 

IF(J] ,fcO,JH2(j) FS(J1.I)=FS(J1,I)/WS0 
IF(J] .EQ.joHj FS ( J1 * I ) =FM 1/WRO 
IF (HflT.LE.FSTOICH) GO TO 820 
010 CONTINUE 

FM?=F { JA, if ) 

JF < Jp,HE. JOH) FM?=Ft1?-Fs (JOH, I)*WTO 
IF (J?.NE. JH2) FM?=FM2-FS(JH2,I) 

IF(J?.NE.JH> FM?=FM2-FS(JH,I) 

IF ( J?.NE, JH20) FH2=Fm2-FS ( JH20? I ) *«G0 
F S(J?.I)=FM2 

IF { J?.EQ. JH20) FS ( J2, 1 ) =FM2/W00 
IF ( J?.LQ. JOH) FS ( J2* I ) =FH2/WT0 
IF (HaT*LE.FSTOICH) Go TO 830 
«20 CONTINUE 

itfr=iter*i 

1 STAR ( I ) = t ter 

IF ( ITFH.lt, IHFGIN) GO TO 689 
C 1 *»««UNDFR-KELaXATION OK MASS FRACTIONS i 

DO 6ns L=1»NS 

HHo FS(L,T)sRPC*KS(LfI)*RPD*STORE(L) 

089 rONTjNUE 

DO 3?5 L=1»NS 

325 RFM (L) = (FS (L. I) -STORE (L) ) / (FS <L» I > +1.E-40) 

C **«LIMITS ON MASS FRACTIONS 
DO 330 Lai , NS 
FS(L, I)sAhIN1(Fs (L»I) *1,) 

339 FS (L, T ) sAhAXI (FS (L, I) *0, ) 

DO 3?6 l.= l,NS 
8? A RT0RE(E)=FS(L,I) 

RMAXsO. 

00 3?Y l.al. NS 

lF(Fs(L,l) ,LT, 1.1-06) GO TO 327 
RMAXbAMaXI (HMAX , A8S (REM ( L ) ) ) 

3?7 CONTINUE 
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C *»*CONVpRGENcE O.'ITtRlON 

IF (ITER.GT.IT MAX. OR »A0S (RMAXJ ... ! .CC) GO TO 208 
GO Tn U5 
?1>8 CONTINUE 

DO 3 AS Jsl.NS 
DO 3 A S I c 1 . NP 1 

345 FS (J. I)=AwAXI (FS(J.I) »0.M 
c ««#tfst FOR POOR CONVERGENCE 
DO 82? f Irlt/voi 

8??J IF(IsTAR(t) .GT.ITMAX) WRITE (6*8828) 1 ST AR ( I ) , I , I STEP 

882* FORMAT ( l0X,« POOR CONVERGENCE <«.I2m ITERATIONS) AT N0DE**I3,« 

1 AND STEP NO.o.ja) 

RETURN 

CHAPTER a 

C*hmhh>LOCA|. TEMPERATURES 
entry eamcla 

DO 9m3 i = 1»NP1 
9M3 PTEMp(l)=F (UTE»I) 

DO AOS 1=1 ,NPI 
SMC=0. 

DO 4pA J = 1 .NS 

rNTH ( I ) = E N T n ( I ) -FS ( J. I ) »HO ( J) 

40 6 SMC = SMC + pi ( J, I ) «CPMN ( J, I I 
405 F ( JTF » I > =pNTH { I ) /SMC 

C #»«UNDEP-REl'aXATION of temperatures in fraction zone 
ILM=iLUCm- 3 
ILP=ILUCP*3 
DO ]A9 1=1. NPi 

IF ( .MOT . ( T .GE.ILV.AND. I.LE. ILP) ) GO TO 149 

F(JTF,I)=.5*(F<JTE t I)fPTEMP<I> > 

149 CONTINUE 
C **#FORM/.T STATEMENTS 
moo format ( 7 f io. 4 ) 

2000 format ( 7 f io.l) 

Return 

end 

COMPILER SPACE 
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SOBRcW* iMr UMTpiiT 

it it it it it it it it it it if it it it it it it it it it it it it it it it it It it it it it it it it it it it it it it it it it it it it it it it it it it it it it 

DIMENSION FLUX (10) , A T I T L 1, (13) ,ATlTL2(l3) , AVRRI..E (6*10) * A3YMRL (7) 
r 0 MM 0 N/GENRAL/ACC.CSALFA,l)PDX»CX.ENTH (4 0) , E ( 1 0 . 4y ) , FS ( 1 0 , AO ) , GAD3 , 
) (iAMM.GAM (/, 0) * I , IFIN* INOE ( 1 0 ) , I ND I ( 1 0 ) , I STEP , IUTRAP * I TEST , K£X , K I N , 
?KRAO,N,IJEo,'IP) ,NM1,OM<40) ,0M(>(40) ,f)()M(40) * K AST (10) , PE I , PS I E , PS 1 1 , 
3P (40) ,R (4r ) *R"o (40) ,RJf < l 0) ,RJ1 (10) « PME ,RMI , SO <40 ) , SU ( 40 ) , WM <1 0 ) , 
4XD,Xn.Y (40) » Yl» IF (40) .YE,YI,RJTE (10) ,RJTI (10) 
rOMMnM/CRflD/RFLO'fl.RFXDtKASE 

COMHOW/Cjc/JUi JK, JO, JHS* JA, JG? JTE, JLE, JUV, JB 
nOMMpH/CF JH2 , JO? , JOH , JH20 , JH * JO * JN2 

C(»MMON/CSTAK/HO (7) ,CP OAR (7,60) ,HT(7,60) ,RC (5,60) , ISTAR (40) ,KA CM 
f:nMMr.N/CTl'»HO/AK , ALMG,CMU.CMU1N,CE1 » CE.2 * CGI ,CG2,CWIT,DUI)Y (40) , 

1 DUD Yen (40 ) .Emul (4(1) , F M U T (4(1 ) , IPO , PR ( 1 0) ,PRT (id) ,RTw (40) 
fOMMCM/OlJTPT/ I AX * I £ NO » Y IK 1 * YOUT , YHA , DYhA , DDYHA » RDUCT * KP , KS » KT , 
J.kASEN ,( o*XsTAT (DO) , XPPOE (20) ,XPLCT (20) , PRESS 
COMMON/TPl.or8/XTAXlS,XTPLOT (40) »YTAXES(10) * YTPLOT (10,40) , 
lYTMAX(lO) . YTSY mR ( 1 0 ) ,OUT<40) ,IPROF 
r0MM0 f J/CPR0P/J02, IN2, Ih20 , OFAC ♦ NR , NS , wCO, WRO , WSO , WTO, GASCON* GAMr.A 

C « it it it it it it it it it it it it it it it itit it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it 

c itit it initial values of interest 
ENTRY OUTP1 

FST0JCH=2.«WM ( JH?) /WM ( J02) 

FSTOiCH^FcTOICH/ (l.+FSTOICH) 

RE AD (5, 1153) (ATITL1 (K) ,K=1,13) 

RE AO ( 5 » 1 1 53 ) (ATITL2(K) ,K = 1,13) - 

READ "(5, 1154) ( (AVRRLE (K.L) »K=1*6) ,|. = 1,10> 

REY = ( HO ( 1 ) «F ( Jij, 1 ) <>2.«Y <NP1) /EPUL( 1 ) 

VM I X e 0 • 

DP 1016 jaltNS 

1016 V M I X = V M I X * E” S ( J ., 1 ) / to M ( J ) 

AMACpaF ( Jii, 1 ) /SORT (GaMmA*GASCON*F ( JTE, 1 ) »VMIX) 

PRESS l=PRrSS 
1 1 1 N = F ( DU , ] ) 

UF IN = F(JU’.1 l-F(JUiNPl) 

0 I N = FS ( JH? , 1 ) 

TIN = F ( JTE". 1 ) 

C iMtitnuTPi'T OF ALPHANUMERIC DATA 

WRITE (6, ) | o2 1 > (ATlfLl (K) ,K = 1 , )3) , (ATITL2 <K) ,K = 1, 13) 

W R I T p ( 6 , i o 2 ii ) (L, ( A VR&LE ( K , L ) »K=1 ,6) ,L=1 ,NEQ) 

V/PITF (6, ln23) 

WRITf (6,ln?4) CPU, CE1 ,CE2, CGI ,CG2 ,AK,aLMG, CHIT, IPD 
WRITE (6, ln3^» ( 0 M ( I ) ,1 = 1, N PI) 

WRITE (6, 1 o33> NR, NS 

READ (5, 1154 > ( (AVRGLE (K,L) , K= 1 , 6 ) , L= 1 , NR ) 

WRITE <6, 1^34) 

WRITF (6, in 22) (|_, (AVR8 lE(K,L) , K = 1 , 6 ) , L= 1 , NR ) 

READ (5, 1154) ( (aVRBLE (K,L) ,K = 1 ,6) ,L = 1 ,NS) 

WRITE (6, j n35) 

WRITE (6, In 22) (L, (AVPRLE ( K , L ) , K= 1 , 6 ) , L = 1 , NS ) 

WRITE ( 6 , ln36) 

WRITE (6, 1 n 25) 

READ (5*1 153) (ASYM8L (K) ,K=1 ,7) 

DO 1 0 1 7 J=l,N(.n 

1017 WRITE (6, ] n2<>) J ,PR ( J ) ,PRT ( J ) , K AST ( J ) 

WRITE (6, In 37) rfy,amach,yin,yowt,rdwct,pre:ssi ,uin,cin,tin 

RE TURN 

r it {.' it it it it it it it it it it it it it it il it it it it it it it it it it it it it it It it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it 
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c ***S t aT ION VftHlAf 5LFS 

ENTRY oujpp 

XUD = yit/Y0llT 

c <hh*mfan square fluctuations (of h elf.mend 
jF(Io?.EO. 0 -AMD.IH 2 O.EQ.n) GO TO 165 
c <u»#locatf position of flame front 

C ASSUMING ^BATTLEMENT* variation of species concentration 

DO 1 A 1 1 = 1* NP 1 

RAT=FST OjrH* ( 1 , -FS ( JN2 , 1 ) ) 

PARA=SFH2-1 .t-10 

SFH2s(F ( Ja ♦ I > -SORT <F < JG, I ) ) ) /RAT 
J F ( Sf'HE.LT . 1 « ) go TO 1 62 
161 CONTINUE 

16? FACIbY(I-1)*(Y(I)-y(I-1) ) « < 1 .-PARA) / (SFh2-PARA) 

no if, 3 i=i,npi 

RAT«FSTOlrH*(l ,-FS(JN2,I ) ) 

PARA=sFH2-l .E-lo 

SFH2s (F ( Ja » I ) +SQRT (F ( JG', I ) > ) /RAT 
IF (SpH2.LT. 1 . ) GO TO 164 

163 CONTINUE 

164 FACE = Y(I-1)*<Y(I)-y(I-1))<‘(1.-PARA)/(SFH2-PARM 

165 CONTINUE 
PRDRp=PRESSl-P ( 1 ) 

ULlHrs (F ( JU» 1 ) -F ( JU*NPl ) )/UEIN 
CLINE=FS ( JH2, 1 ) /(CIN+1.E-30) 

TLINe=F (JTE*1) /TIN 

WRITE <6, In06) I STEP , , XUO , OX . YHa , 0 YUA , DDYHA , YOX 

URITf (6, lnO'H K IN .KEX » YIN , Y (NP1 ) , POE ,ACC > CPU » CEE 
WRITE f 6. in 05) P S 1 1 » PS I E » P E I * RM I *RME * FAC I * FACE »F ACM 
WRITE (6, 1 oO^) F (JU,1 ) »F (JA,1) ,F ( JTE * 1 ) » ULINE , CL InE , TLINE , PROWP 
C «##CH£CK for flux conservation 
DO 1 n?7 J=liNF0 

FLUX | J) = (F { J»?) *0M (2) *F (J.N)# (1 .-OM(N) ) ><»2. 

DO ln?7 I*2*NM1 

10 27 FLUX { J) =FLU-<(J) + <F ( J, I \ *F ( J, 1*1 ) ) *0M0 ( I ) 

00 ln?H J=1 »NI-q 

IF (KaST (J) ) 1 0 ? 9 * l 0 2 9 * 1 0 2 6 

1029 F LUX'(J)=FLUX(J)-,25»( (F (J.3) -F< J»2J ) * (GM (3) -QM (2) )♦ 

1 (F ( J,NMl ) -F ( J,N) ) « (ON (N) -OM(NMl ) ) ) 

10 2A FLUX '< J) s. 5«ELUX { J) «PFUPSI I«F ( J, 1 ) -PSXE«F ( J.NPl ) 

WRITE ( 6 , In30) (J.ELUX (J) , J = 1,NEQ) 

IF (I<ASE.Eo. 2) WRTTf (6*1155) RFLOW , RDUCT , REXD , P ( 1) , DPDX 
IF (KeX. NE. l) RETURN 
OF =P JF t JHc ) /H ( NP 1 ) 

TAUE=RJE ( JIJ ) /R ( NP 1 ) 

WRITE ( 6 * 1 n 1 0 > JNDE(JHS) ,F(UTE,NP1) .QE.TAUE 
return 

c #«« •**<!<> <MHH> « (HHHHUt# •» QQHQ # tt « tt « ft « » ft ftftft ft 0 ft ft ft ft ft ft ft ft ft ft ft ftftft 

ENTRY 0UTP3 

c «<HtpRoF i t.E Variables 



C/ INFORMATION ", TEMPORARY) on new output routine. 
c/ smoroutine prof i l assigns values for plot, ano also writes 
C/ PROFILES. IT IS CALLED H Y ... 

C,/ r ALL PRUFtL ( JPROP. TITLE. FIRST. ADO, D IV, final, XPLOT. SYMBOL) 

c/ whfpf jpruebj refers to The fij.n array 

C/ JPRfiF .gT • 0 MEANS USE F(UPROF.I) ARRAY AND WRITE PROFILE. 

C./ JPRCiF . ) T » 0 MEANS USE OUT H)* ARRAY AND WRITE PROFILE. 
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C/ JPRoF=0 MEANS NO ACTION' UNLESS KPLOT.LT.O 
C M.IT. OUT (I) ARRAY IS OVERWRITTEN BY PROFIL. 

C/ TITLE IS TH'r Name OF THE PROFILE, PLOTTED BY SYMBOL. 

C/ FIPS? ^NU FINAL ARE FIRST ANO LAST VARIABLES WHITTEN IN PROFILE 

C/ Ann AND OIV modify PROFILES from 1=2 TO I=N 

C/ rPLOT.OT.O means assign YTPLOT (KPLOT, I ) ARRAY. 

C/ «PL OT-O MEANS NO PLOT ASSIGNMENTS. 

C/ VPLOT.LT. 0 MEANS ASSIGN THE Y ( I ) ARRAY TO THE XTPLOT(I) ARRAY. 

C/ SYMBOL TS THE CHARACTER USEO IN PLOT. 



C «n»*ZERO YTPLoT ARRAYS 

IF (ISTEP .GT. 1) GO To *01 
DO 40? 1=1, NP1 
40? nUTm=0. 

DO 40? Ks 1 1 1 0 

403 CALL PROFtL( 0,1H«,0.,0.,1,,0., K,lh*) 

A 0 1 rONTlNUE 
C *<k>DISTaNCES (Y) 

CALL PROFfL( i* * 6 H Y / Y ( N P • R (1) » 0 . , Y { NP 1 ) , Y ( NP 1 ) , - 1 , 0 ) 

C **»AXIAl. VELOCITIES (JU) 

DIV = F (JU,1 )-F ( JlJ.NPl) +1.E-30 

CALL PRRFtL ( JU.GHU Vr.LO , F ( JU , 1 ) , -F ( JU , NP1 ) , D I V , F ( JU * NP1 ) , 1 , 1 HU ) 

C ###STaGI"ATION ENTHALPIES (JHS) 

()JV = F (JUS’. 1 )“F (JHS,NPl) 

CALL PROFrL ( JHS.6HH ST aG . F ( J riS , I ) ,-F(JhS»NPl) , 0 I V , F ( JHS , NP1 ) ,0,0) 

C ###TliRU».'LENT KINETIC ENERGIES (JK) 

0IV=(F(JU’. I)-F(Jii,NPI)+1.E-30)«»2 

CALL PROp ,L ( J*< , AHK .E . T|) ,F ( JK , 1 ) , 0 . , Dl V, F ( JK ,NP 1 ) , 5»1HK) 

C <n»»0)SSIPATlcN RATES (JD) 

CALL PROFrL < JD.fiHO 0 I SS , F ( JO , 1 ) , 0 . , 1 . , F ( JO , NP 1 ) , 6 , 1HO ) 

C *#«OISSIPATIcN LENGTH SCALES (JLE) 

00 l?nl l=p,N 

1201 F ( JLF, I ) =<sQHT ( F ( JK , I ) ) / ( RTW ( I ) + 1 . E-30 ) 

f t jlf , l ) . 

F(JLF,NPl)aO, 

DIV=Y (NP1 ) 

CALL PROFfL(JLE,GHLENGTH,F{JLE,I) , 0 , , 0 1 V , F { JLE , NP I ) , 10,1HL) 
c *<M>^TURf’"LENT VISCOSITIES^ (EMUT) 

'••'RITf (O.ftifi) (EMOT(I) ,1 = 1, NP 1 ) 

C **«DENSITIES 

WRITE (^»6r*U) (RHO ( I ) , 1 = 1 *NP1 ) 

C «<h»REyN0!.US STRESS CORRELATION (JUV) 

OIV=(F(JU'.1)-F(JU,NP1)+1.E-30)<*«2 

CALL PROFjL ( JUV.ftHijV 1 , F ( JljV , 1 ) ♦ 0 . , 01 V , F ( JUV , NP I ) , 2 , 1HS ) 
c <nn>sPECiES mass fractions (ja> 

DIV=F (OA, i ) -F (JA,NP1) ♦l.E-JO 

CALL PROFrL <JA,7HF(JA,I) , F ( JA , 1 ) , -F ( JA , NP 1 ) , D I V , F < JA , NP 1 ) , 3 , 1HA ) 

concentration fluctuations cof h element) 

0IV = F (JA, l ) -F <JA,NP1 ) +1 .E-30 

CALL PR0FiL(JG, 7HF(JP,I) ,F(JG,1) ,-F(JG,NPl) ,01 V ,F ( JG,NP1 ) ,9«1HG) 

C #«#ABSO(JlTF TEMPERATURES (JTE) 

niV = F(JTE'. 1 ) -F (JTE. NPl ) +1. .E-30 

CALL PKOFtL (JTe .AHTf.MPiF < JTE, 1 ) , -F ( JT.E.NP 1 ) , 0 1 V , F ( JT£ , NP 1 ) , M , 1 n V ) 

C «MH>OIITPilT OF INDIVIDUAL SPECIES CONCENTRATIONS 
DO ftp* Jal.NS 

A?5 WRITE (G,t,?6) ASYI'BL ( J) , <FS ( J, I ) , 1 = 1 ,NI>1 ) 

WRITE (G,6*>7) (ISTAH(I) ,I = I ,NPl) 
c »#»PLOTTlMG nF Pr>t)FJ LES 
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605 IF(IsTEP.FO.O) RFTURN 

IF ( IsTLP.rO. IAX.OM. ISTEP.EQ. IEND) GO TO 435 
IF (Xrn.LT.XPLOT (KT) ) RETURN 

kt=kt*1 

435 wRiTF(bflo9l) XudSTEP 

CALL. PLOT*? (XTPLOT,35.NPl,XTAXIS*YTPUOT,YTMAX,iO,10,YTAXES,YTSY Mb) 

C tt ************* tt **************** «HHM> * *« **£•»■»*«» tt# ***»**»«»#•,>»*>« <HH). <>****</ ••! 
c ***FORMAT STATEMENTS 

6?4 FORMAT <1H , A6 , 1 P l 1 E 1 1 • 3/ < 7 X , 1 Q 1 1 F 1 1 . 3 >> 

6?7 format < ih ,6hitek ,nin/<7x, mil) ) 

6?ft FORMAT ( 1 H * 6HEMIJT * 1 P 11 E 1 1 . 3/ ( 7X , 1 P 11 E 1 1 . 3 ) ) 

6?9 format ( in .ghrho »ipiieii.3/<7x*ipiieii.3) > 

100a FORMAT (//1HV,*I$TEP=**I5. 1X,*XU=*,1PE11.3,1X,*XUD=*> 1PE1 1 .3* IX, *0/ 
1=*. lpFll .3,1 X,*yHA=*. 1 PE 11.3, 1X*«DYHA=*,1PE11.3, 1 X , *DDYHAs* , IRE 11 . 
23 * 1 X , *YDX=* » 1PE 1 1 .3 ) 

10 07 FORMA T ( 1 x >K I N = * , I 2 * ix , «KEX=* , I 2 , 1 X , * Y I N=* , 1 PE 1 1 . 3 , 1 X , *Y ( NP i ) = * , 1 F 
1F11.3- 1X,*PUE=*,1P£11.3,1X,#ACC=*, IPEll . 3 , 1 X , *CMU=* , 1 PE 1 1 . 3 , 1 X , *CE 
??=tt,1PEli f 3) 

10 0ft FORMAT < 1 X.*PSII=*,1PF 11.3, IX, *PSIE=*,1PE11.3,1X,*P£I=*,1PF1 1.3, IX. 
]*RMIs»*1PfU . 3; 1X,*RME=*, 1PE11.3, 1X,*FACI=*, 1PE11.3, lX.*FACb=*i 1PE 
211 .3, lX,*FACMs« f 1PE11.3) 

10 09 FORMA T <) X - . *U1 =*, 1 P£l 1.3, 1 X, *M1 =<S) PEI 1. 3, IX, *T1=*, IPEll. 3,1 X,*Ul/i 
1 i=#, J.PE11 .3, 1X,*M1/MI=*. 1PE11 .3, 1X,*T1/TI=*, IPEll .3t IX , *PRDHP=* . 1- 
2F11 .3) 

1010 FORMAT (1X>INOe ( JHS) = <>, 12 » l X , *TW aLL=* » 1 PE 1 1 . 3 , IX , *HT . TRANSF . =* . 1P£ 
111 ,3,1.***sHEAR STRESS=»,1PE11.3) 

1021 FORMAT (1H1 ,4X,13A6/5X,13A6///5X»» ORDER OF DEPENDENT VARIABLE 

IS IS*/) 

102? FORMAT (6X', I2»*. *,6a6) 

1 o?3 FORMAT <//lH0,*— CONSTANTS IN TURBULENCE MODEL* ) 

1 0 ?4 FORMAT </4X,»CMU=*,F8. 3,3 X,*CE1=*, FB. 3,3 X,*CE2 = *,F8. 3, 3X , *CGl ^ . 
.)3.3X,*CG2=*»FB.3,3 X,«Ak=*,FG.3,3X,*ALMG=*,F8.3,3X»*CRIT=*,FH.3*3x, 

2*IPDj:*, 12) 

1025 FORMAT (/1H0»*-- PRANDTL SCHMIDT NUMBERS ARE*/ ) 

10 24 FORMAT (4X',*U=*, I2,3X,*PH(J) =*,F9.3,3X,*PRT (J) =*,F9.3,3X, *K AST (J) =* 
1.12) 

loan FORMaT(5h flux,, 7(1H{,I2,2H)=,1Pe10.3,2X) ) 

1032 FORMAT ( iHn omega DISTRIBUTION — «//(4X,lPllEll.3) ) 

1033 FORMAT (/1H0«*- NUMBER OF REACTIONS =*,I3,10X,*- NUMBER OF SPECIES 
1=*,I3) 

1034 FORMAT (/1H0.5X,* — REACTIONS ARE*) 

1 035 FORMAT (/1H0,5X,#— CHEMICAL SPECIES ARE*) 

1034 FORMAT <////) 

1037 FORMAT </lHO»*— INITIAL VALUES — */5X,*REY =*, IPEll »3/5X,*MACH =* 

1 dPEil ,3/5X,*YlN =*, IPEll. 3/5 X,*VOUT =*, 1PE1 1 .3/SX,*R0UCT=* , lPt 1 1 

?.3/5X, tt PRrSSs*, IPEll, 3/5X»*U =* * 1 PE 1 1 . 3/5X ♦ *C =* , 1PE1 1 . 3 /Sa , 

3*T a * , 1 PE 1 1 . 3 ) 

1091 FDRMaT(4h1XU=,1PE10.3,8H, ISTEP=,I6> 

1153 Forma T ( 1 3a 6 ) 

1 1 54 Format <6Aa ) 

1155 format ( i x>rflou = *, IPEll .3 , ix,«rduct=*, IPEll .3, ix,*kexu=*. ipki i .a. 
] 1 X » *P°ESS=* * 1 PE 1 2 » 4 » 1 X * *DPOX=* * 1 PE 1 1 • 3 ) 

C********************** **************************************** ********* 
RFTUpN 

FMD 

COMPTI.FP SPACE 
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S' •BrtfllT tl'Jp t’HOFlL ( JPROF * T I TL E , FIRST * ADD , I) I V , F iNAL * KPLOT , SYmSO L) 
COMMoN/GpNR AL/ACC , CS ALF A , QPDX ; DX , ENTM (40) * F (10*40) , FS ( 1 0 , 40 ) ; GAM 3 , 
lRAMM,r,AM(«0) * I * IF IN* INDE ( 10) ,INDl<10) ♦ I STEP , I UTRAP * I TEST * KEX ,klM, 
?kRaD,M*NEo»'^M ,NMJ. ,OM (AO) t OMD ( 40 ) ,ROM (40) * K AST ( 1 0 ) . PE I * PS I E * PS I I * 
3P (40) • R ( 4 n ) * RHO (AO) tRJE(lO) , R J I (10) * RME * RMI * SO ( 4 0 ) , SU ( 40 ) * WM ( 1 0 ) * 
4X0, Xu, V ( 4 n ) * VO IF (40) , YE*Yl ,RJTE ( 1 0 ) , p JT I { 1 0 ) 

COMMPN/CJC/JU. JK. JD* JHS* JA, JG * JTF,JLE, JUV, JB 
COMMON/CF<;/JHp. jnp, JOH, JH/?0, JH, JO* JN2 

COMMON/TPLOTM/XTAXIS,XTPLOT(40) *YTAXES(10> «YTPLOT (10*40) * 
lYTMAx ( 10) . YTSYOB ( 10) *0UT (40) , IPHOF 

C 

C«»<hmk>sUQHo'JTINf To WRITE PROFILES and ASSIGN TRANSVERSE PLOT ARRAYS 



IF(Kp!.OT) 500,700*700 
700 IF(JPROF) 200*200,100 

inn DO lnl 1=1 *NPI 

101 OUT ( I ) B F ( JPROF , I ) . 

200 YMAX=1 • E-3 0 

DO 20 1 1=1 *NP1 
YMAXbAMAXI (YMAX,0UT (I ) ) . 

201 OUT ( T ) = (OUT < I) ♦ ADO) /l)IV 

IF ( JpROF.NE.O. AND. IPROF.FO. 1 ) WRITE (6,900) TITLE, FIRST, (OUT ( I) » 1 = 2 
3 , N ) * FINAL 

IF (JPROF. NF..O. AND. IPROF.E^.2) WRITE (6,900) TITLE, (F (JPROF, I) ,1=1 ,N 
JP1 ) 

TF(KPl.OT) 500.600,400 

400 DO 40 l I=l,NPl 

401 YTPLOT(KPLOI,I)=OUT(I) 

YTAXfS (KPLOT) sTITLE 
YTSYl'H (KpLOT ) s.SYMBOL 
YTmAx (KPlpT) =ymax 
return 

500 on 5n 1 1 = 1, NP1 
: sol XTPLoT (I) = (Y ( [) +ADD)/0IV 
XTAXIS=TITLE 

write (t»,9oO) TITLE, FIRST, (XTPLOT (I) ,1=2, N >, FINAL 
aoo Return 

900 FORMaTUH ,A6,1P11E11.3/(7X,1P11E11.3) ) 

END 


COMPILER SPACE 



SUBROUTINE PLOTS (X.IDIm, I MAX, XAXIS, Y,YMAX, JDIM, JMAX, YAXES, SYMBOL) 



C SUBROUTINF F ot? PLOTTING J CURVES OF Y ( J , I ) AGAINST X(i) 



C 

C X A NO V ARE ASSUMED TO BE' IN ANY RANGE EXCEPT THAT NEGATIVE VALUES 

C ARE plotted AS ZERO 

C X AND Y ARE SCALED to THE RANGE 0. TO I. BY DIVISION BY THE MAXIMA 

C WHICH are also printed 

C I D I M IS THE VARIABLE DIMENSION FOR X, 

c tmax is the number or x values 

C XAXIS STORES THE NAME OF THE X-AXIS 

C JD I M jS THE VARIABLE OIMENSION F()R Y. 

C JMAX IS THE NUMBER OF CURVES TO BE PLOTTEO, (UP TO 10). 

C THE ARRAY YAXES(J) STORES THE NAME OF THE CURVES. 

C THE ARRAY SYMBOL (J) STORES THE SINGLE CHARACTERS USED FOR PLOTTING 

C .. . 

DIMENSION X ( IDJM) ,Y (JDlM.IDIM) ,YMAX (JDIM) »YAXES (JDIM) , SYMBOL ( JDIM ) 
l.A ( 10J ) 

DATA qOT , rROSS, BLANK/1 H, , 1 H + , 1H / 

C*«hmh> SCALING X ARRAY TO THE RANGE 0 TO SO 
XUA*el.E-30 
DO 1 1=1,1 MAX 

1 IF (X'(T) .GT.XMAX) XMAX = X(I) 

PO 2 1=1, TMAX 
X ( I ) = X ( I ) /XMAX4»50. 

? TF (X (T ) .LT.O.) X ( I ) =0. 

SCALING V ARRAY to the RANGE 0 TO 100 
DO 3 J= 1 , JMAX 
AYMAX=l.E-30 
DO 4 1=1 , TMAX 

4 AYMAXsAMAXl UYM aX,Y(J»I)) 

DO 3 1=1, TMAX 

3 Y ( J, T ) =AMaX 1 { ( Y { J, I ) /AYMAX)#lO0. ,0. ) 

€*>»»»« IDENTIFYING THE VARIOUS CURVES TO BE PLOTTED 
WRITE (6, ln3) XAXIS, XMAX 
WRITE ( 6 , 1 o 0 ) (YAXES(I) ,1=1, JMAX) 

WRITE (G,ln6> (SYMBOL (I) ,1 = 1, JMAX) 

WRITE (6, ln2) (YMAX(I) ,1=1, UMAX) 

DO 5 1=1,11 

5 A(i)=0.1#FLOAT(I-l) 

WRITE (6, lol ) (A ( I ) , 1 = 1 , 1 1 ) 

MA I N LOOP. EACH PASS PRODUCES AN X-CONSTANT LINE 
DO 4o 1 = 1 ’,51 

IF(I.E«.1.0R.I.E0.51) GO TO 32 
GO TO 33 

ALLOCATE . OR ♦ AS MARKER ON THE Y-AXIS 
3? DO 30 K = l’, 101 

30 A (K)=DOT 

DO 3} K=u, 101,10 

31 A (K ) =CROS«: 

C«*hmm*aLLOCaTE . OR ♦ MARK ON THE X-AXlS, ALSO THE APPROPRIATE X VALUt. 

33 A(l) e nOT 
A ( 1 0 ] )=DOT 
, K= I- 1 
4A K=K-5 

IF (K) 4fl , 7 , 4B 
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47 A(i)=C w OS<; 

A (ion a CknSS 
4a XL=0,n2«FLOAT< I - 1 ) 

check if ANY Y ( X ( I ) ) value lies ON THIS X-CONSTANT line 
IF YES G n TO 41. OTHERWISE GO TO 42 
00 43 K=i'. IHAX 

IF(IFI*CX#K)*1.5)-I> 43.41.43 
C*<hhh> L OCATE y ( K(I) ) 

41 DO 44 J = l'.JMAX 
NY = Y ( J.K) *1 .5 

a (NY) s symbol < j> 

44 CONTINUE 
GO To 4 2 
43 CONTINUE 

PRINT X-rONSTANT LINE 
4? WRlTr (6, ln5) XL, (A (K) ,K=1 ,101 ) .XL 
putting blanks into x-constant line 

00 49 K=r. 101 

49 4(K)sRLANK 
■ 49 CONTINUE 

oo 5n I = r.u 

50 A(I)s.l*FLOAT (1-1) 

WRITF (9.104) (A (J) ,1 = 1,11) 

RETURN 

100 FORMAT ( 1 1H Y-AXES ARK » 5X , 1 0 (1 X » A 1 0 ) ) 

101 FORMAT (1H0.2X.11F 10.1) 

105 FORMAT (1SH MAXIMUM VALUES . 1P1 <»£ 1 1 . 3 ) 

103 FORMAT (11H0X-AXIS IS .A0.17H .MAXIMUM VALUE s.1PE10.3) 

104 FORMA T<3X'. I IF 10. 1/1 HI) 

105 Format <2H X » F6 . 2. 3X ♦ 1 01 A1 »F6.2) 

10k F0RMaT(7h SYMBOL, 11X, 10(1X, A10) ) 

FND 

compilpr SPACE 
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SUBWOIJTINf YIWT (FRaC.YFRAC, JJ) 

COMMoN/GeNRAL/ACC,CSaLFA,0PDX,DX,ENTH(40) ,F (10,40) ,FS{ 1 0,40) * GAM3 « 
1GAMN,GAH<;;o) *I«IFIN«lKDE(ia) ,INRI (10) ,ISTEP,IUTRaP,ITEST,KEX7m.N, 
?KRAD,N*IMfO»NPl ,NM1 ,0m (40 5 ,OMD<40) ♦ BOM (40) ,KAST(l0) , PE I , PS I E » PS 1 1 , 
3P (40) (40) (40) ,RJE (1 0) »RJI (10) , RME »RMI , SO ( 40 ) , SU (40 ) , WM ( 1 0 ) , 

4 X |.) « X I ; » Y (4ft) ,YI)IF(40) . YF » Yl *RJTE <1 0 ) ,RJTI (10) 

COMMON/Cjc/JU, JK, JD, JHS, JA, JG» JTE, JLE 1 , JUVtvlB 
C0MM0N/CFS/JH2, J02, JOH, JH20, JH, JO, JN2 



C<n»#*oiNTEf.'POLATlON SUBROUTINE 

C-t 

KEFOIF=FRaC»ab$ ( f (JJ,NP1)-F<JJ,1) ) 

IF (WFFOIF.NE.O. ) GO TO 10 
y-FRACcPRAroY (NPl ) 

WRITE ( 4 , 1 ) YFPAC 

1 FORMAT 17h lSTEP=,I4*4H XU= , 1PE1 0 . 3 , 52(1 IN SUBROUTINE YINT , HEF0iF= 
10, SO YFRaC=ERAC<*Y(NP1)=,E10,3) 

WRITE (0,2) FRAC,JJ,NP1»F(JJ*NP1> ,F(JJ,1) 

? Format (6h frac=, ipeio.3,4h jj=,I3,5h npi?*,I4, 

1 11H F(JJ,NP1)=,E10.3,9H F(JJ,1)=*E10.3) 

return 

lo IF (FPAC.GT. .5) GO TO 30 
DO 20 I =2 i, N 

DIF=APS(F (JJ,I)-F (JJ, 1) ) 

IF (DIF-REFOIE) 20,21,21 
21 T=ABs(F(jj,I)-F(JJ,I-l) ) 

TF(T) 22,33,22 
?3 T= (DIF-REFDIF)/T 
23 YFRAC = Y(I)-T*(Y(I)-Y(I-D) 

RETURN 
20 CONTINUE 
30 On 4n IDA<;H = 2,N 
I=N-,?-1DacH 

OIF=ARS(Ff JJ, I ) -F ( JJ, 1 ) ) 

IF (DJF-REFOIF) 41 , 41,40 
41 T = AHS(F(JJ,1)-F(JJ,I*1') ) 

IF (T) 42,43*42 
4? T= (UIF-REFDIF) /T 
43 YFRAC=Y(I)-T«(Y(I*1)-Y(I) ) 

RETURN 
40 CONTINUE 
RFTUPN 
END 

compiler SPACE 
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